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ABSTRACT 

A  local  geoid  model  to  predict  the  geoid  heights  in  the  vicinity  of  Monterey 
Bay,  California,  was  developed  to  use  Global  Positioning  System  (GPS)  differential 
positions  and  known  Mean  Sea  Level  (MSL)  with  the  method  of  collocation.  The 
local  geoid  models  were  based  on  Rapp's  360  degree  x  360  order  global  geoid 
model  determined  from  gravity  measurements.  Control  data  were  adjusted  by  least 
squares  to  solve  for  the  parameters  in  the  local  geoid  model.  Also  studied  were 
factors  that  affected  the  GPS-measured  ellipsoid  height  differences.  These  included 
(1)  comparing  GPS  differencing  solutions,  (2)  standard  error  of  GPS  observations, 
(3)  corrections  for  surface  meteorological  values,  and  (4)  observation  durations  for 
GPS. 

The  data  used  in  this  research  were  taken  from  GPS  measurements  on  the 
campus  of  the  Naval  Postgraduate  School  (NPS),  an  area  about  100  m  x  630  m  and 
in  an  area  approximately  15  km  x  33  km  near  Monterey,  California.  The  time 
period  was  from  February  5,  1988,  to  May  12,  1988. 

The  accuracy  of  the  predicted  geoid  heights  is  ±  2  cm  if  a  six-parameter  model 
is  used  for  the  large  area,  and  ±  2  to  10  mm  if  a  five-parameter  model  is  used  for  the 
NPS  campus. 
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I.   INTRODUCTION 

The  Global  Positioning  System  (GPS)  is  able  to  establish  precise  relative 
positions  in  the  World  Geodetic  System  of  1984  (WGS  84).  A  Trimble  GPS 
receiver,  which  has  the  capability  of  measuring  carrier  phase,  was  used  to 
determine  the  vector  base  line  in  space.  The  components  of  the  base  line  are 
expressed  in  terms  of  cartesian  coordinate  differences  (AX,  AY,  AZ)  [Remondi, 
1984].  These  vector  base  lines  can  be  converted  to  distances,  azimuths  and  the 
ellipsoid  height  differences,  Ah,  relative  to  the  WGS  84  Ellipsoid. 

The  results  of  several  tests  and  operations  have  clearly  shown  that  GPS 
survey  methods  can  replace  conventional  horizontal  survey  methods. 
Comparable  accuracies  have  also  been  achieved  for  GPS-derived  ellipsoid  height 
differences,  Ah.  The  problem  of  converting  these  ellipsoid  height  differences, 
Ah,  to  orthometric  height  differences,  AH,  remains  to  be  resolved.  For  example, 
in  engineering  surveying  applications  WGS  84  coordinates  must  be  transformed 
to  the  North  American  Datum  of  1983  (NAD  83)  system.  The  GPS  obtains 
ellipsoid  heights,  rather  than  the  orthometric  heights;  the  geoid  height,  N,  must 
be  calculated  to  obtain  the  latter.  One  of  the  problems  in  this  transformation  is 
the  accurate  determination  of  the  local  geoid  height,  N. 

For  GPS  survey  applications,  a  geoid  model  should  provide  geoid  heights 
with  an  accuracy  commensurate  with  that  of  the  ellipsoid  height,  h,  so  the 
accuracy  of  the  derived  orthometric  heights  is  not  reduced.  In  the  future 
differential  positioning  will  be  widely  used  in  GPS  surveying,  and,  therefore, 
only  the  geoid  height  differences  between  stations  will  be  required. 

Geoid  height  computation  techniques  include  the  following: 

(1)  Geoid  height  differences  in  the  U.S.  can  be  determined  from  gravity 
data  and  the  Stokes'  integral  method,  or  from  astrogravimetric  data  and  least 
squares  collocation  methods.  These  methods  lead  to  uncertainties  that  are 
typically  1  to  10  cm  for  distances  of  less  than  20  km  and  5  to  20  cm  for  distances 
between  20  to  50  km  [Zilkoski,  1988] 

(2)  A  comparison  of  a  data  set  from  GPS  with  gravimetrically  determined 
geoid  heights  using  least  squares  collocation  techniques  shows  discrepancies 
between  the  two  data  sets  of  about  ±  2  cm  for  a  maximum  intersection  distance  of 
approximately  50  km  [Denker  and  Wenzel,  1987]. 
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(3)  Mean  gravity  anomalies,  deflections  of  the  vertical  and  a  geopotential 
model  calculated  to  degree  and  order  180  have  been  used  to  determine  geoid 
heights  in  the  area  bounded  by  (34°  <  <p  <  42°,  18°  <  X  <  28°)  [Tziavos,  1987]. 
The  method  used  was  that  of  least  squares  collocation.  By  using  empirical 
covariance  functions  for  the  data,  suitable  weighting  functions  for  the  different 
sources  of  observations,  and  the  optimum  cap  radius  around  each  point  of 
elevation,  an  accuracy  better  than  ±  0.60  m  was  obtained  for  geoid  heights. 

The  main  objective  of  this  thesis  is  to  use  a  model  that  predicts  N  in  a  region 
near  Monterey,  California,  from  the  GPS  differential  positions  and  known  mean 
sea  level  (MSL)  using  the  method  of  collocation.  Also  studied  were  factors  that 
affected  the  ellipsoid  height  differences,  Ah,  obtained  from  GPS  measurements, 
such  as  comparing  GPS  differencing  solutions,  standard  error  of  GPS 
observations,  corrections  for  surface  meteorological  value  ,  and  observation 
durations  for  GPS  measurements.  The  data  used  in  this  research  were  taken 
from  GPS  measurements  on  the  campus  of  the  Naval  Postgraduate  School  (NPS) 
an  area  about  100  m  x  630  m  and  in  an  area  approximately  15  km  x  33  km  near 
Monterey,  California.  The  time  period  was  from  February  5,  1988,  to  May  12, 
1988. 

The  accuracy  of  the  predicted  geoid  height  is  ±  2  cm  if  a  six-parameter 
model  is  used  for  the  large  area,  and  ±  2  to  10  mm  if  a  five-parameter  model  is 
used  for  the  NPS  campus. 


II.   GEOID  HEIGHT 


A.    THE  GEOID 

Surveyors  and  engineers,  in  most  cases,  are  interested  in  the  orthometric 
height,  H,  as  measured  above  the  reference  surface  of  the  geoid.  The  ocean  is 
considered  to  be  freely  moving,  homogeneous  and  only  subject  to  the  force  of 
gravity.  When  a  state  of  equilibrium  is  achieved,  the  surface  of  this  idealized 
ocean  assumes  a  level  surface  of  the  gravity  field.  It  may  be  regarded  as  also 
extending  under  the  continents.  This  level  surface  is  called  the  geoid.  If  the 
potential,  W,  is  given  as  a  function  of  the  coordinates  r,  §  and  X,  then  the  geoid  is 
given  by  [Moritz,  1984] 

W  =  W(r,(j)^)  =  W0 

The  geoid  is  a  closed  and  continuous  level  surface  which  extends  partially 
inside  the  solid  body  of  the  earth.  The  direction  of  the  gravity  vector  at  any 
point  (plumb  line  or  vertical)  is  normal  to  the  geoid.  The  curvature  of  the  geoid 
displays  discontinuities  at  abrupt  density  variations.  Consequently,  the  geoid  is 
not  an  analytic  surface,  and  therefore  not  a  practical  reference  surface  for 
position  determinations.  The  geoid  however,  is  well  suited  as  a  reference  surface 
for  potential  or  height  differences,  which  are  obtained  by  precise  levelling  in 
combination  with  gravity  measurements. 

To  establish  the  geoid  as  a  reference  surface  for  heights,  the  ocean  water 
level  is  recorded  and  averaged  over  long  intervals  (>  1  year)  using  tide  gauges. 
The  MSL  thus  obtained  represents  an  approximation  to  the  geoid.  The  National 
Geodetic  Vertical  Datum  of  1929  (NGVD  29)  was  derived  for  land  surveys  from 
a  general  adjustment  of  the  first  order  levelling  net  of  both  United  States  of 
America  and  Canada.  In  the  adjustment  MSL  was  observed  at  twenty-one  tide 
stations  in  the  United  States  and  five  in  Canada.  The  geoid  established  by  this 
method  may  deviate  by  ±  1  to  ±  2  m  from  a  level  surface  due  to  periodic, 
nonperiodic  and  secular  variations  [Torge,  1980]. 


B.    THE  WGS  84  ELLIPSOID 

The  development  of  WGS  84  [DMA,  1987]  was  initiated  by  the  United  States 

Department  of  Defense  for  navigation  and  weapon  systems.     The  Defense 

Mapping  Agency  (DMA)  developed  WGS  84  as  a  replacement  for  WGS  72.  The 

defining  parameters  and  reference  frame  orientation  of  the  WGS  84  Ellipsoid 

and  the  WGS  84  Ellipsoid  Gravity  Formula  are  those  of  the  internationally 

sanctioned  Geodetic  Reference  System  of  1980  (GRS  80).    Accordingly,  a 

geocentric  equipotential  ellipsoid  is  defined  by  the  semimajor  axis  (a),  the 
flattening  (f),  the  equatorial  gravity  (ya),  and  the  angular  velocity  (co).  The  WGS 

84  Ellipsoid  used  the  values: 


a 

= 

6378137    m 

f 

= 

1/298.257223563 

Ya 

= 

987.03267714    gals 

co      =      7.2921 15xl0"5   rad  /  s 

The  reference  system  for  GPS  is  WGS  84.     The  precise  geocentric 
coordinates  obtained  from  GPS  receivers  are  in  WGS  84. 


C.    THE    ORTHOMETRIC    HEIGHT,    H,    THE    ELLIPSOID 

HEIGHT,  h  AND  THE  GEOID  HEIGHT,  N 

In  geodetic  applications  three  different  surfaces  or  earth  figures  are 
normally  involved.  First  is  the  earth's  actual  topography;  second,  the  geometric 
surface,  or  ellipsoid  and;  third,  the  equipotential  surface,  the  geoid.  The 
relationship  between  the  earth's  topography,  the  ellipsoid  and  the  geoid  in  a 
section  through  the  earth's  surface  is  shown  in  Figure  1 . 


NORMAL  TO  GEOID 


NORMAL  TO  ELLIPSOID 


Figure  1.  Relationship  between  the  earth's  surface,  the  geoid  and 
the   ellipsoid. 

Two  features  in  the  figure  are  of  particular  interest  : 

1.  The  deflection  of  the  vertical,  0,  defined  by  Pizzetti  as  the  angle  at  the 

geoid  between  the  direction  of  the  plumb  line  (normal  to  geoid)  and 
the  normal  to  the  ellipsoid  through  the  point,  P,  on  the  geoid  [Torge, 
1980]. 

2.  The  vertical  separation,  N,  between  the  geoid  and  the  ellipsoid. 


The  deflections  of  the  vertical  and  the  geoid  heights,  N,  depend  on  the 
ellipsoidal  coordinates  and,  hence,  on  the  parameters  of  the  reference  ellipsoid 
and  its  position  with  respect  to  the  earth.  If  they  are  referred  to  the 
geocentrically  situated  mean  earth  ellipsoid,  then  they  are  referred  to  as  absolute 
quantities;  otherwise,  they  are  relative  quantities.  The  absolute  deflections  of  the 
vertical  in  flat  terrain  and  the  highlands  assume  values  between  one  and  ten 
seconds  of  arc;  in  mountainous  areas,  they  vary  between  one-half  and  one  minute 
of  arc.  As  a  result  of  density  variations  within  the  earth,  the  geopotential 
surfaces,  including  the  geoid,  have  irregular  shapes.  Absolute  geoid  heights, 
however,  rarely  exceed  100  m  [Torge,  1980]. 

The  orthometric  heights,  H,  shown  in  Figure  2  are  referred  to  an 
equipotential  surface,  the  geoid.  The  orthometric  height  of  a  point  on  the  earth's 
surface  is  the  distance  from  the  reference  surface  to  the  point,  measured  along 
the  plumb  line  normal  to  the  geoid.    The  ellipsoid  height,  h,  of  a  point  is  the 
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distance  from  the  reference  ellipsoid  to  the  point,  measured  along  the  line  which 
is  normal  to  the  ellipsoid.  For  purposes  of  simplicity,  the  orthometric  height,  H, 
and  the  ellipsoid  height,  h,  are  shown  to  be  along  a  common  vertical.  In  most 
cases,  this  would  cause  a  very  small  error  that  is  considered  insignificant 
compared  to  present  uncertainties  of  the  geoid  height  N  estimates.  The  geoid 
height,  N,  is  defined: 

N  =  h  -  H 
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Figure  2.     Relationship   between  the   geoid  height,  N,  the  ellipsoid 
height,  h  and  the  orthometric  height,  H. 


The  orthometric  height  has  greater  physical  meaning  than  the  geometrical 
ellipsoid  height.  The  orthometric  height  has  traditionally  been  determined  by  the 
technique  of  levelling  in  which  increments  of  height  are  obtained  from  the 
intersection  of  the  line  of  sight  of  a  level  instrument,  tangential  to  the 
geopotential  surface  and  passing  through  the  level  axis,  with  two  graduated  rods. 
Accurate  orthometric  height  information  is  needed  for  precise  engineering 
operations  such  as  the  construction  of  dams,  pipelines,  and  tunnels. 

Geoid  heights  have  been  accurately  determined  for  some  major  geodetic 
datums  such  as  the  North  American,  European  and  Australian.  These  datums  are 


well  supplied  with  astrogeodetic  deflections  and  have  fair  gravity  coverage.   The 
standard  error  of  relative  geoid  heights  in  these  areas  is  about  2  or  3  m. 

D.    THE  GEOID  HEIGHT  FROM  GRAVITY  MEASUREMENT 

Both  the  U.S.  National  Geodetic  Survey  (NGS)  of  National  Oceanic  and 
Atmospheric  Administration  Service  (NOAA)  and  Trimble  Navigation  versions 
of  Rapp's  360  degree  x  360  order  model  (OSU86F)  were  available  to  me.  Rapp 
computed  two  potential  coefficient  fields  that  are  complete  to  degree  and  order 
360  [Rapp,  1986].  One  field  (OSU86E)  excludes  geophysically  predicted 
anomalies,  while  Rapp's  other  model  (OSU86F)  includes  such  anomalies.  These 
fields  were  computed  using  a  set  of  30-minute  mean  gravity  anomalies  derived 
from  satellite  altimetry  in  the  ocean  areas  and  on  land  from  standard 
measurements. 

Gravity  anomalies  can  be  observed  and  then  used  to  compute  the  geometric 
deviation  of  the  geoid  from  the  ellipsoid.  The  expression  for  the  gravitational 
potential  is  written  in  the  following  form  [Rapp,  1986]: 

kM  °°  Tali     l    (  \ 

V(r,(b,X)  =  1  +  I    -      X      C,  cosm>.  +  5,    sinmk      P,  (sin<b) 

Y      L  1=2  LrJ     mt0V       I™  lm  J       l™  . 

where 

r,  (J),  X  :  geocentric  coordinate 

kM  :  geocentric  gravitational  constant 

a  :  equatorial  radius  of  the  reference  ellipsoid 

C,  ,  S,  :  fully  normalized  potential  coefficients 

lm       lm 

P       :  fully  normalized  Legendre  function  of  degree  1  and  m 

y       :  normal  gravity 
The  potential  at  a  point,  U,  is  the  scalar  sum  of  V  and  centrifugal  force  potential 
V  [Ewing,  1976]: 

U(r,<j),?i)   =  V(r,<ja)  +   V'(r,(j),?i) 

The  difference  between  observed  gravity  potential  W  and  the  computed  normal 
gravity  potential  U  is  denoted  by  T  [Moritz,  1984],  so  that 


W(r,fcA.)  =  U(r,<j),?i)  +  T(M>,X.) 


compared  the  geoid  defined  by  the  potential  W  is  given  by 


W(r,<|>,A.)  =  W, 


A  reference  ellipsoid  with  the  same  potential,  W    =  U  ,  is  given  by 
U(r,4>,A.)  =  W0 

A  point  P  on  the  geoid  is  projected  onto  the  point  Q  on  the  ellipsoid  along  the 
normal  to  the  ellipsoid  PQ  =  N.  PQ  is  the  distance  between  geoid  and  ellipsoid  at 
the  point  and  is  called  the  geoid  height  (Figure  3). 


GEOID 


W  =  W 


o 


ELLIPSOID 

u  =  w0 


Figure  3.    Geoid  and  ellipsoid  (Moritz  [1984,  Fig.  2-12]). 


The  disturbing  potential  T  can  be  written  as  [Rapp,  1986]: 


T(r,fcX)  = 


kM  ~ 

~v"X 

y  1=2 


1  I    I   CtYf^) 

m=0  a=o 
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where 


C£r   {^,01  =  0,  and  S    ,o=l) 


Y^((f),^)  =  {P   (sin(j))cosmA.,  a  =  0,  and    P   (sin<|))sinm^,  a  =  1 } 

Since  the  N  is  relatively  small  compared  to  the  ellipsoid  geocentric  radius  R', 
then  [Ewing,  1976] 

T  =  yN 

and 

T 

.\   N  =  - 
Y 

where  y  is  the  normal  gravity  at  (<j),  X). 

The  gravity  anomaly,  Ag,  in  the  Molodensky  surface  free-air  anomaly  sense 
is  given  by  [Moritz,  1984]: 

gp   • ■   gp   "  Yq 

The  boundary  condition  that  relates  Ag  and  T  is  [Moritz,  1984] 

a  9T      1  3T_ 

As  =  -ah+F  3hT 

where  h  is  the  distance  along  the  plumb  line  direction.  Neglecting  deflections  of 
the  vertical  we  have  [Rapp,  1986] 


Ag=^-I(I-D 


x  £(  ££,- eft -32- 32  )*&♦.*) 


Y       1=2  Lrj   m=0  a=0 

where    5^  ( i  =  h,  y)  are  ellipsoidal  corrections. 

The  Stokes  function,  S(\\f),  can  then  be  used  to  solve  the  geoid  heights  above  the 

geodetic  ellipsoid  [Ewing,  1976]. 


—Jo  "Jo*  Ag(\|/,a)S(y)sin\(«i\(Kia 


Ns=  Ji 


where 

Y        :  the  mean  value  of  normal  gravity 

m 

V|/  :  the  angular  distance  between  the  point  where  N  is  being 
determined  and  the  area  where  the  effect  of  Ag  is  being 
considered 

a        :  the  azimuth  from  the  affected  point  to  that  causing  the  effect 

Ag      :  the  gravity  anomaly 

S(y)   :  the  Stokes  function 

The  Stokes  function  is  given  by  [Ewing,  1976] 

S(\j/)  =  esc  J-  +  1  -  6  sin  5r-  5  cosy  -  3  cosy  ln(sin^+  sin2  £) 

The  gravitational  potential  using  degree  360  and  order  360  has  been  tested 
through  comparison  of  Doppler  station  geoid  heights  with  geoid  heights  from 
Rapp's  versions  of  geopotential  models.  The  agreement  between  the  two  geoid 
height  measurements  is  approximately  ±  1.6  m  [Rapp,  1986]. 
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III.   GEOID  HEICxHT  FROM  GPS 


A.    THE  GLOBAL  POSITIONING  SYSTEM 

The  Global  Positioning  System  (GPS)  calls  for  a  precise  navigation  system 
divided  into  three  segments:  space  segment,  control  segment  and  user  equipment. 
The  space  segment  will  consist  of  three  orbital  planes  of  satellites  at  inclinations 
of  120°  in  circular  orbits  at  altitudes  of  20,000  km.  Each  plane  will  eventually 
contain  six  to  eight  satellites  to  give  the  three  dimensions  of  position,  velocity  and 
precise  time  24  hours  a  day  anywhere  in  the  world.  The  control  segment  consists 
of  the  ground  stations  necessary  to  track  the  satellites,  monitor  the  system 
operation,  and  periodically  provide  corrections  to  the  navigation  and  time 
signals.  Each  satellite  broadcasts  signals  containing  information  on  its  position. 
The  GPS  satellite  transmits  signals  at  two  L-band  frequencies  (1227  and  1575 
MHZ)  to  permit  corrections  for  ionospheric  corrections.  The  signals  are 
modulated  with  two  codes:  P,  which  provides  for  precise  measurement,  and  C/A, 
which  permits  easy  lock-on  to  the  desired  signal.  The  user  segment  consists  of 
the  equipment  necessary  to  convert  the  satellite  navigation  message  into  useful 
navigation  information. 

The  navigation  message  contains  the  data  that  the  user's  receiver  requires  to 
perform  the  operations  and  computations  for  successful  navigation  with  GPS. 
The  data  includes:  (1)  information  on  the  status  of  the  Space  Vehicle  (SV);  (2) 
time  synchronization  information  for  the  conversion  of  the  C/A  to  P  code;  (3) 
parameters  for  computing  the  clock  correction  and  the  ephemeris  of  the  SV;  and 
(4)  corrections  for  delays  in  the  propagation  of  the  signal  through  the 
atmosphere.  In  addition  the  data  contain  almanac  information  to  define  the 
approximate  ephemeris  and  to  give  the  status  of  all  the  other  SV  information 
which  is  required  for  use  in  signal  acquisition  [Milliken,  1980].  Ranges  to  the 
satellites  are  determined  by  signal  transit  times  multiplied  by  the  speed  of  light 
(299,792,458  m/sec).  The  transmitted  message  contains  ephemeris  parameters 
that  enable  the  user  to  calculate  the  position  of  each  satellite  at  the  time  of  the 
transmission  of  the  signal. 
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B.  CARRIER  PHASE  MEASUREMENTS 

GPS  measurements  can  be  made  using  the  pseudo-range  and  the  carrier 
phase.  The  pseudo-range  is  essentially  a  measurement  of  distance  contaminated 
by  clock  error.  When  four  satellites  are  observed  simultaneously,  the  three 
dimensional  position  of  the  ground  receiver  can  be  determined  along  with  the 
receiver  clock  offset  at  a  single  epoch.  The  accuracy  of  pseudo-ranges  is  affected 
by  multipath  effects,  which  depend  on  the  antenna  design,  and  its  height  above  the 
ground. 

Carrier  phase  measurements  are  more  precise  than  pseudo-range  and  are 
not  as  vulnerable  to  multipath  effects.  They  can  be  used  to  compute  the  precise 
base  line  components  AX,  AY,  AZ  between  two  receivers.  Phase  measurements 
are  made  by  beating  the  received  carrier  with  the  signal  from  a  local  oscillator 
internal  to  the  GPS  receiver.  The  slant  range  from  a  GPS  receiver  to  a  satellite 
can  be  modelled  in  terms  of  time.  It  takes  the  signal  time  to  travel  between  the 
satellite  and  the  receiver  or  the  equivalent  number  of  cycles.  The  range  of  the 
cycles  will  consist  of  an  integer  and  fractional  number  of  cycles.  When  a 
receiver  locks  onto  the  carrier  signal,  it  can  immediately  measure  the  fractional 
part  and  begin  counting  subsequent  integer  cycles,  but  it  can  not  measure  or 
account  for  the  initial  integer  number  of  cycles  that  preceded  the  initial  fractional 
part.  The  initial  integer  ambiguity  which  biases  the  subsequent  measurements  is 
called  the  initial  integer  ambiguity  bias. 

GPS  uses  a  one-way  carrier  beat  phase.  The  GPS  satellite  and  receiver  are 
controlled  by  separate  clocks.  The  satellite  clock  generates  the  signal,  and  the 
receiver  clock  detects  when  the  signal  arrives.  An  error  in  the  synchronization 
of  the  clocks  of  1  microsecond  creates  an  error  in  range  of  300  m. 

C.  ONE-WAY       CARRIER       PHASE       MEASUREMENT 

DIFFERENCING 

A  single  difference,  SD(j,i),  is  formed  by  differencing  carrier  beat  phase 
observables  from  two  receivers  1  and  2  at  the  same  observation  epochs  i  of  same 
satellite  j.  The  equation  is  given  by  [Remondi,  1984]: 

SD(j,i)  =  S(2,j,i)-S(l,j,i) 
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where  S(k,j,i)  is  the  raw,  unpreprocessed,  fractional  phase  plus  the  count  made  at 
epoch  i  by  receiver  k  for  satellite  j.  The  main  advantage  of  the  single  difference 
is  that  it  reduces  or  eliminates  satellite  orbital  and  clock  errors,  because  they  are 
common  to  both  receivers.  Its  disadvantage  is  that  one  can  not  exploit  the  integer 
nature  of  integer  ambiguities.  Thus,  for  short  base  lines,  the  ultimate  in  accuracy 
may  not  be  achievable  [Remondi,1985]. 

A  double  difference,  DD(j,k,i),  is  formed  by  differencing  single  differences 
between  a  reference  satellite  j  and  another  satellite  k  at  the  same  epoch  i.  The 
equation  is  given  by  [Remondi,  1984]: 

DD(j,k,i)  =  SD(k,i)  -  SD(j,i) 

The  advantage  of  the  double  differences  is  that  the  receiver  clock  dependent 
terms  are  eliminated  because  the  differences  for  each  epoch  are  correlated.  The 
significance  of  the  removal  of  those  terms  is  to  reduce  from  nanoseconds  to 
microseconds  the  timing  accuracy  required  to  achieve  one  cycle  accuracy.  For 
short  base  lines,  the  integer  ambiguities  can  be  isolated,  since  the  contribution 
made  by  the  clock  drift  is  reduced  [Remondi,  1985].  The  Trimble  4000SX 
receiver  achieves  sub-microsecond  accuracy  by  using  the  C/A  code  timing 
information  [Ashjaee,  1985]. 

A  triple  difference,  TD(j,k,i),  is  formed  by  differencing  the  double 
differences  for  the  same  satellite  pair  at  some  integer  of  succeeding  epochs  i+1 
[Remondi,  1984]. 

TDG,k,i)  =  DDQJci+1)  -  DDG,k,i) 

The  advantage  of  the  triple  difference  is  that  it  eliminates  all  the  time  independent 
terms,  namely  the  initial  integer  ambiguities,  and  becomes  insensitive  to  the 
initial  ambiguities  and  any  cycle  slips  when  the  receiver  loses  lock.  The 
disadvantage  of  the  triple  difference  is,  another  level  of  correlation,  loss  of 
resolution  and  a  greater  number  of  observations.  Triple  differences  are  already 
correlated  with  respect  to  satellite  because  of  the  underlying  double  differences 
and  are  further  correlated  with  respect  to  time  because  consecutive  triple 
observations  will  have  common  DD(j,k,i+l)  terms. 
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For  short  base  lines  local  area  integer  ambiguities  can  easily  be  resolved 
because  unmodelled  errors  are  highly  correlated  between  the  two  antenna  sites  and 
are  mostly  eliminated  by  differencing.  Algorithms  can  take  advantage  of  the 
integer  nature  of  the  initial  ambiguities  and  solve  for  them  [Remondi,  1984]. 

D.    GEOID  HEIGHT  FROM  GPS 

Surface  fitting  techniques  can  be  used  with  GPS-derived  geoid  heights.  GPS 
stations  are  likely  to  be  close  together,  of  the  order  of  a  few  tens  of  km,  and  the  local 
geoid  can  be  estimated  directly  if  levelling  data  is  available  in  the  area.  This 
together  with  it's  relative  simplicity  makes  the  method  practical  for  correcting  GPS 
heights. 

It  is  possible  to  use  the  two  sets  of  elevations  (that  is,  levelled  and  GPS- 
derived)  to  define  two  distinct  planes.  The  published  levelled  elevations  are 
referred  to  the  geoid  and  the  GPS  elevations  are  referred  to  an  ellipsoidal  surface. 
If  both  elevations  are  made  equal  at  one  bench  mark  (control  station),  then  in 
general,  the  other  bench  marks  will  have  two  elevation  values.  After  several 
different  models  were  studied  two  mathematical  surface  models  were  chosen  in  this 
study.  The  five-parameter  model  found  suitable  for  use  on  the  NPS  campus  is 

N0(X,Y,Z)  =  h0  +  Ah  -  H  +  aAY  +  bAX2  +  cAY2  +  dAXAY 

and  the  six-parameter  model  selected  for  a  larger  area  near  Monterey,  California,  is 
given  by 

N0(X,Y,Z)  =  h0  +  Ah  -  H  +  a'AX  +  b'AZ  +  cAX2  +  d'AY2  +  eAXAY 


where 


N0(X,Y,Z) :  the  global  geoid  height,  obtained  by  gravity  measurements 

h0     :  the  ellipsoid  height  of  the  reference  point,  including  a  constant 

correction  to  N0(X,Y,Z)  at  the  reference  station 
Ah    :  the  ellipsoid  height  difference  with  respect  to  the  reference  station 
H     :  the  published  or  levelled  elevations  at  each  control  station 
a,  b,  c,  d,  a',  b',  c',  d',  e' :  the  coefficients  to  be  determined 
AX,  AY,  AZ  :  the  coordinate  differences  in  WGS  84 
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These  models  can  be  solved  to  determine  the  east-west  and  north-south  tilts 
which  are  absorbed  by  two  rotations,  one  around  the  north  axis  (AY)  and  the 
other  around  the  east  axis  (AX)  in  the  horizon  system.  Their  separation  is 
absorbed  by  the  scale  correction  [Zilkoski,  1988].  The  local  geoid  heights,  N, 
can  be  found  from  the  equation 

N  =  N0(X,Y,Z)  +  AN 

where  AN  is  the  variation  of  geoid  height  in  local  area.  The  equation  is  given  by 
AN  =  H  +  N0(X,Y,Z)  -  Ah 


or 


AN  =  h0  +  aA  Y  +  bAX2  +  cAY2  +  dAXA  Y 


AN  =  h0+  a'AX  +  b'AZ  +  c'AX2  +  d'AY2  +  e'AXAY 


To  solve  for  h0,  a,  b,  c,  d,  a',  b',  c',  d'  and  e',  the  global  geoid  heights 
N0(X,Y,Z)  can  be  obtained  from  the  Geoid.exe  program  described  in  Chapter 
IV.  The  geoid  model  can  be  rearranged  to  give 


or 


H  =  h0  +  Ah  -  N0(X,Y,Z)  +  aAY  +  bAX2  +  cAY2  +  dAXAY 

H  =  h0  +  Ah  -  N0(X,Y,Z)  +  a'AX  +  b'AZ  +  c'AX2  +  d'AY2  +  e'AXAY 


Then  h0  a,  b,  c,  d,  a',  b',  c',  d'  and  e'  can  be  solved  by  the  least  squares  method. 
The  geoid  height,  N,  found  by  this  method  appears  to  be  adequate  for  areas  up  to 
50  km  x  50  km  where  the  geoid  is  smooth  [King,  1985]. 

E.    METHOD  OF  COLLOCATION 

The  method  of  collocation  was  derived  from  least  squares  interpolation. 
This  method  was  used  to  predict  the  geoid  height,  N,  for  a  local  area.  The  geoid 
model  is  given  by 

N  =  N0(X,Y,Z)  +  AN  +  S  +  n 
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where 

S  :  the  signal 

n  :  the  noise 

N0(X,Y,Z)  +  AN  :  the  system  function  from  the  surface  models  at  control 

points  where  both  h  and  H  are  known 
Then,  for  a  given  N0(X,Y,Z)  the  method  of  collocation  can  be  used  to  determine  S 
at  these  control  points  and  to  predict  S  at  other  points  in  the  local  area.  At  any  point 
in  the  local  area  the  value  of  h  can  be  determined  by  using  differential  GPS 
measurements  between  a  known  point  and  any  other  point  by  the  equation 

h  =  h0  +  Ah 

The  value  of  H  at  any  point  is  given  by 

H  =  h    -  N0(X,Y,Z)  -  AN  -  S  -  n 

The  value  of  H  at  any  point  in  the  local  area,  which  depends  on  the  value  of  h  at 
the  control  points  and  the  differential  GPS  measurements,  can  be  predicted  with  an 
accuracy  of  ±  n.  The  general  form  of  the  observation  equation  in  the  method  of 
collocation  is  [Jeyapalan,  1977]: 

x  =  A»X  +  S  +  n  +  0«S 

q      q         p 


where 


x  :  the  vector  of  the  observation  (x  =  Ah  -  N0(X,Y,Z)  -  H) 

A  :  a  given  rectangular  coefficient 

X  :  the  vector  of  the  systematic  parameters  (h0,  a,  b,  c,  d,  or  h0,  a',  b',  c\ 

d\  e') 

S  :  a  signal  vector  at  q  observation  points 

n  :  a  vector  of  measuring  errors,  noise  at  q  points 

S  :  a  signal  vector  at  p  unknown  stations 

•  :  indicates  matrix  multiplication 

0  :  the  null  matrix 


If 


Z   =  S  +n 

q        q      q 
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Then 

x  =  A^X  +  Z  +  0«S 

q  q  p 

then 

X  =  (ATCq"1A)-1  ATCq"1x 
and 

S   =  C    C  -1  (x  -  A  X) 

p         pq    q 

where  the  variance-covariance  matrix  C  of  the  Z  and  S   vectors  is 

q  p 

(C      C     \ 
c  -       p       pq 
u  "    C       C 

V       pq  q  J 

The  essence  of  this  method  is  that  by  some  means  a  covariance  matrix  can 
be  assigned  to  the  signal.  For  noise  it  will  be  possible  to  assign  a  diagonal  weight 
matrix. 

S  are  the  values  of  the  signal  at  the  interpolated  stations.  Suppose  there  are 
q  observations  (and  values  of  S  ),  p  interpolated  values  of  S  and  m  model 
parameters.  The  covariance  matrices  are  the  following  [Bomford,  1980]: 

(i)  C  ,  the  expected  covariance  between  the  observed  x's  for  all  pairs  of 

the  q  observations.  It  is  a  q  x  q  matrix. 

(ii)  C   ,  the  expected  covariance  between  the  signals  for  all  pairs  of  the  q 

observations.   It  is  a  q  x  q  matrix. 

(iii)  Cnq,  the  expected  value  of  n   at  each  station.   It  is  a  q  x  q  diagonal 
matrix. 

(iv)  C      the  expected  covariance  between  all  pairs  of  mixed  observed 

and  interpolated  signals.  It  is  a  p  x  q  matrix. 

(v)  C  ,  the  expected  covariance  between  the  signals  at  pairs  of 

interpolated  stations.  It  is  a  p  x  p  matrix. 

The  variance  of  the  noise,  C  can  be  estimated  in  the  usual  way  according  to 
the  circumstances  at  each  station,  different  types  of  instrument,  etc.  The  noise  at 
different  points  is  independent;  hence 
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c 

nq 


/^an0 0  ^ 

0  o 


lo  O....0J 


where  a„  is  the  standard  error  of  the  noise 
by 


The  expected  covariance  between  the  observation  x's,  C  can  be  obtained 


C    =  C     +  C 

q  nq  sq 

because  the  signal  is  small. 

The  C    can  be  be  computed  from  a  simple  function  whose  parameters  can 

be  determined  by  using  empirical  data.  The  covariance  function  is  positive  and 
definite.  These  two  characteristics  are  found  in  many  functions.  Functions 
commonly  used  for  covariance  may  be  constant,  sinusoidal,  Gaussian, 
exponential,  exponential  cosine,  and  exponential  sine  and  cosine.  In  this  thesis  the 
sinusoidal  function  was  used.  The  function  is  given  by 

C     =  B  sin(kr) 

sq 

C(r)  =  C     +  B  sin(kr) 

nq 

where  Cnq  is  the  standard  deviation  of  the  control  stations  and  B,  k  are 
coefficients  to  be  determined. 

The  parameters  can  be  determined  by  least  squares,  and  residuals  at  each 
point  can  be  computed.  From  the  residual,  the  covariance  between  points  at  a 
distance  (r)  can  be  computed  by 

L*     0    1 

cm  =  - — r 

v  y         n  -   1 

where  V0  is  the  residual  at  the  center,  Vr  is  the  residual  at  a  point  which  is  at  a 
distance  r  from  the  center  and  n  is  the  number  of  points.  There  are  several 
methods  of  determining  C(r)  of  which  the  concentric  circle  approach  is  the 
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simplest.    The  coefficients  B  and  k  can  then  be  computed  from  the  computed 
covariance. 

The  expected  covariance  between  the  signals  at  pairs  of  interpolated  stations 
C   can  be  obtained  by  substituting  the  distance  between  the  interpolated  stations 

and  control  stations  into  the  covariance  function. 

The  expected  covariance  between  all  points  of  mixed  observed  and 
interpolated  signals  C     can  be  obtained  from  the  covariance  function  for  each 

distance  from  the  interpolated  stations  to  the  control  stations. 

In  this  thesis  two  cases  were  studied  in  solving  for  the  geoid  height  model. 

Case  1.  Cq  =  I    and  0^  =  0 
Case  2.  Cq  *  I    and  C^  *  0 

I  start  with  assuming 

C    =  P  *  =  I 


C     =  0 

pq 


where  P  is  the  weight  matrix  and  I  is  the  unit  matrix. 
Then 


and 


XQ=(ATPA)4  ATPx 


S   =  0»P(x-AXJ  =  0 
p  o 


Using  X  to  compute  residuals,  v,  and  then  estimating  C  ,  C     and  P,  it  is  then 

q       pq 

found  that 


X=(ATC"1A)-1  ATC  *x 
1  a  a 


and 


q 
-i 


S  =C    C"'(x-AX) 

p  pq      q  V 
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IV.   DATA   COLLECTION  AND  INSTRUMENTATION 


A.    PLANNING 

Nine  temporary  bench  marks  were  established  on  the  NPS  campus  (Figure  4). 
Bench  marks  GH7  and  GH8,  designated  as  check  marks,  were  established  in  the 
center  of  the  levelling  loop.  H  of  bench  marks  was  obtained  by  differential 
levelling.  Ah  was  measured  with  a  pair  of  GPS  receivers. 


A  TREE 

N 

S  antenna 

\ 

_      Main  Gate 


A  GH6 


A  GH7 


▲  GH5 


10  m 


Figure  4.    The  temporary  bench  marks  on    the  NPS  campus. 
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Eight  permanent  bench  marks  were  recovered  in  the  study  area  (Figure  5) 
[Vertical  Control  Data,  1961].  Bench  mark  GWM  27,  designated  as  a  check  mark, 
was  roughly  in  the  center  of  the  survey  area  for  the  permanent  bench  marks.  H  was 
obtained  from  the  published  elevations.  Ah  was  also  measured  with  a  pair  of  GPS 
receivers. 


Watsonville 


121°  45' W 

+ 


121°  32' W 

+    36°53'N 


N 


▲  J  697 


▲  D  697 

Salinas 


GWM  27  A  p  gl2 
▲  S  812 


+     36°  38'  M 


*  NPS 


Monterey 


10  km 


Figure  5.    The  permanent  bench  marks  in  the  Monterey  Bay  area. 
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B.    DATA  COLLECTION 

1.   Obtaining  H 

a.    Precise  Levelling 

Precise  levelling  was  run  on  the  NPS  campus  on  February  19,  1988, 
and  March  24,  1988.  The  elevation  of  station  TREE  was  assumed  to  be  zero 
above  the  MSL,  so  the  orthometric  height  differences,  AH,  for  each  temporary 
bench  mark  could  be  observed.  The  differences  of  the  relative  elevations,  AH, 
for  the  forward  and  the  backward  sights  on  the  NPS  campus  are  listed  in  Table  1. 


Table  1.  AH  FOR  THE  NPS  CAMPUS 

Bench  mark 

AH(m) 

From 

To 

Forward 

Backward 

Mean 

TREE 

GH1 

-0.808 

0.808 

-O.808 

GH1 

GH2 

5.388 

-5.384 

5.386 

GH2 

GH3 

3.150 

-3.149 

3.149 

GH3 

GH4 

1.535 

-1.536 

1.535 

GH4 

GH5 

-1.016 

1.016 

-1.016 

GH5 

GH6 

-4.214 

4.215 

-4.214 

GH6 

GH7 

0.489 

-0.489 

0.489 

GH7 

TREE 

-4.031 

4.033 

-4.032 

GH6 

GH2 

0.547 

-0.546 

0.547 

GH7 

GH2 

0.057 

-0.057 

0.057 

GH3 

GH8 

0.716 

-0.718 

0.717 

GH8 

GH4 

0.817 

0.817 

0.817 

To  avoid  an  obstruction  near  bench  mark  D  697,  which  is  10  ft  north 
of  a  12-inch  cedar  tree,  the  temporary  bench  mark  TEMP  1  was  established 
nearby.  Offset  levelling  was  run  on  April  8,  1988.  The  elevation  of  TEMP  1 
was  offset  by  -0.401  m  from  the  elevation  of  D  697. 

F  813  which  was  set  in  the  west  end  of  the  south  abutment  of  a  steel 
bridge  over  the  Salinas  River  was  offset  to  temporary  bench  mark  TEMP  2.  The 
offset  levelling  was  run  on  April  9,  1988.  The  elevation  of  TEMP  2  was  offset 
by  -2.218  m  from  the  elevation  of  F  813. 
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b.   Published  Elevation 

H  values  for  the  larger  of  our  two  study  areas  in  and  near  Monterey, 
California,  areas  were  obtained  from  the  published  elevations  printed  by  U.S. 
Department  of  Commerce  Coast  and  Geodetic  Survey  Washington  D.C.  [1961]. 
The  geodetic  datum  used  in  this  publication  was  the  NGVD  29.  First-order  spirit 
levelling  has  extended  this  datum  over  most  of  the  continent.  Although  first- 
order  lines  may  be  300  km  apart  in  some  western  areas,  most  points  in  the 
country  are  no  more  than  50  km  from  an  estimated  first-order  bench  mark.  A 
readjustment  of  the  this  network  is  the  NAVD  88.  This  new  adjustment  will  be 
made  to  the  geopotential  surface  rather  than  to  sea  level,  and  it  will  place  the 
existing  vertical  data  in  a  form  that  makes  it  most  consistent  and  accessible  to  the 
user.  Changes  to  older  published  elevations  are  not  expected  to  exceed  15 
decimeters  [NASA,  1978].  Table  2  lists  the  orthometric  heights  of  the  permanent 
bench  marks  we  used. 


Table  2. 
STUDY 


H  FOR  PERMANENT  BENCH  MARKS  UESD  IN  THIS 


Bench  mark 

H(m) 

K152 

16.965 

B21 

5.787 

S812 

15.978 

P812 

13.134 

D697 

30.057 

J  697 

85.061 

F813 

7.983 

GWM27 

63.487 

2.  Obtaining  Ah 

a.    Satellite  Observation  Plan 

Due  to  the  positions  of  the  satellite  orbits  during  this  study,  the 
observing  window  was  between  2100  and  0200  hours  Pacific  Standard  Time 
(PST=UTC+8  hours.).  The  time  period  was  between  February  2,  1988,  and  May 
12,  1988. 
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b.  Satellite  selection 

The  same  five  satellites  (SV)  (6,  9,  11,  12  and  13)  were  used  for  all 
observations. 

c.  Position  Dilution  of  Precision  (PDOP) 

The  accuracy  to  which  positions  are  determined  using  GPS  depends 
on  two  factors:  (i)  satellite  configuration  geometry,  and  (ii)  measurement 
accuracy.  GPS  measurement  accuracy  represents  the  combined  effect  of 
ephemeris  uncertainties,  propagation  errors,  clock  and  timing  errors,  and 
receiver  noise. 

The  effect  of  satellite  configuration  geometry  is  expressed  by  the 
dilution  of  precision  (DOP)  factor,  which  is  the  ratio  of  the  positioning  accuracy 
to  the  measurement  accuracy  [Wells,  1987]. 

o  =  DOP  •  o0 

where 

c0  :  is  the  measurement  accuracy  (standard  deviation) 

o    :  is  the  positioning  accuracy  (standard  deviation  in  one  coordinate) 
The  value  of  GDOP  itself  is  a  composite  measure  that  reflects  the 
influence  of  satellite  geometry  on  the  combined  accuracy  of  the  estimation  of 
observation  time  (user  clock  offset)  and  receiver  position  [Milliken,  1986]. 

GDOP2=    PDOP2     +  TDOP2 

TDOP  is  the  Time  Dilution  Of  Precision,  the  error  in  the  clock  of  the  receiver 
bias  multiplied  by  the  velocity  of  light.  The  four  best  satellites  selected  by  the 
receiver  are  those  with  the  lowest  GDOP.  Trimble  recommends  that  the  rapidly 
changing  PDOP  provides  better  geometry  for  phase  differencing  techniques. 
Low  or  constant  PDOP  provides  weaker  solutions.  The  4000SX  receiver  does 
not  record  GDOP,  but  it  does  record  PDOP  every  five  minutes.  PDOP  was 
about  4.5  for  all  GPS  observations. 
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C.    INSTRUMENTATION 

1.  Levelling 

A  Zeiss  model  Ni-2  level  (Ser.  #  82377)  was  used  at  the  temporary 
bench  marks  for  the  third-order,  class  I  levelling.  It  has  a  32-diameter 
magnification,  produces  an  erect  image,  and  has  stadia  constants  of  333  or  100. 
A  Peg  test  was  performed  before  levelling.  The  level  error,  or  c-value  was  also 
checked  before  the  beginning  the  levelling.  The  c-value  was  -0.006  mm/m  which 
was  less  than  +0.05  mm/m,  so  it  was  not  necessary  to  adjust  the  level  [Bodnar, 
1975].  The  level  contains  a  bubble  tube  to  permit  positioning  parallel  to  the 
geoid.  When  properly  set  up  at  a  point,  the  telescope  is  locked  so  that  it  will 
rotate  through  a  360°  arc  in  a  horizontal  plane.  With  the  level  locked  in  position 
readings  are  made  on  two  calibrated  staffs  held  in  upright  positions  ahead  of  and 
behind  the  instrument.  The  difference  between  readings  is  the  difference  in 
elevation  between  points.  Dietzgen  Model  #  6450  metric  rods  were  used  in  the 
levelling.  These  rods  are  graduated  in  centimeters.  The  actual  reading  is 
estimated  to  the  nearest  1  mm.  Rod  levels  were  used  to  indicate  when  were 
vertical. 

2.  £PS 

a.   Trimble  4000SX  Receiver 

A  complete  description  of  the  Trimble  4000SX  receiver  is  given  by 
Trimble  Navigation  [Trimble,  1987a].  NPS  operates  three  Trimble  4000SX  GPS 
receivers.  For  this  study  I  fixed  one  antenna  to  the  roof  of  Building  224  on  the 
NPS  campus,  and  one  was  carried  to  the  field.  The  4000SX  is  capable  of 
observing  the  C/A  code,  integrated  Doppler,  and  carrier  beat  phases  of  up  to  five 
satellites  simultaneously.  Its  ability  to  use  the  C/A  code  allows  the  receiver  to  be 
used  as  a  stand  alone  navigation  system  to  determine  positions  using  Doppler- 
smoothed  pseudoranges  and  velocities  [Ashjaee,  1985].  The  receiver  uses  the 
C/A  code  in  a  time  transfer  mode  to  determine  the  offset  and  drift  of  its  own 
clock  and  thus  provide  accurate  time  tags  for  the  observations  without  the 
requirement  of  an  external  atomic  clock  or  synchronization  with  the  receiver  at 
the  other  end  of  the  baseline.  The  reference  position  (the  geodetic  coordinates  of 
the  antenna)  and  the  practical  options  chosen  must  be  entered  into  the  receiver  via 
the  receiver  key  pad.    The  4000SX  requires  115  V  AC  power  or  12  V  DC 
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power  supply.    115  V  AC  power  was  used  for  the  NPS  campus  measurements. 
12  V  DC,  supplied  by  a  car  battery,  was  used  in  the  field. 

b.  Grid  Personal  Computer  (PC) 

For  precise  relative  positioning  the  4000SX  receiver  transmits  data 
through  an  RS-232  port  to  a  microcomputer  (Grid  PC)  for  storage  on  floppy 
disks  for  post-processing.  The  Grid  PC  uses  either  115  V  AC  or  12  V  DC.  On 
the  NPS  campus  a  115  V  AC  power  supply  was  used,  so  measurements  could  be 
made  for  four  hours.  The  12-volt  battery  in  the  Grid  PC  lasts  about  120 
minutes,  so  measurements  were  taken  for  only   100  minutes  in  the  field. 

c.  Antennas 

Multipath-resistant  Trimble  microstrip  antennas  were  installed  on 
Building  224  on  the  NPS  campus  and  over  the  various  bench  marks.  The  antenna 
heights  from  the  center  of  bench  marks  to  the  edge  of  the  antenna's  ground  plane 
were  measured  before  and  after  GPS  observations.  The  field  antenna  was 
mounted  on  a  tripod  with  a  tribach  and  optical  plummet  for  centering  and 
levelling  the  antenna.  Arrows  on  the  antennas  were  oriented  to  the  north  at  both 
stations  using  a  magnetic  compass. 

d.  Meteorological  Instruments 

A  barometer  and  a  sling  psychrometer  were  used  to  measure 
atmospheric  pressure,  relative  humidity  and  air  temperature  at  each  field  site. 

e.  The  steel  tape 

A  three-meter  steel  tape  was  used  to  measure  the  antenna  height. 
This  tape  can  be  read  to  1  cm.  Readings  are  estimated  to  1  mm. 

D.    SOFTWARE 

A  complete  description  of  the  Trimble- supplied  Trimvec  software  is  to  be 
found  in  Trimble  Navigation  [Trimble,  1987b].  The  software  provides  data 
logging,  baseline  computation  and  datum  transformation  programs.  These 
operate  with  an  IBM  compatible  personal  computers.  A  printer  and  a  hard  disk 
are  used  with  the  planning  and  processing  programs.  The  following  programs 
are  used  on  3-1/2  inch  micro-floppies: 

1.   Data  Logger  on  Disk  1 

In  relative  positioning,  the  receiver  is  controlled  from  the  Grid  PC  by 
version  D  of  Trimble's  Gridlog5.bat  program.  Each  observation  session  was 
initialized  to  log  data  when  a  minimum  of  four  satellites  were  15°  above  the 
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antenna's  horizon.  Five  satellites  were  designated  for  each  observing  session. 
The  observables  and  receiver  clock  parameters  were  logged  to  a  floppy  disk 
every  15  seconds  and  the  C/A  code-determined  antenna  position  and  PDOP  every 
five  minutes.  The  GPS  navigation  message  was  logged  to  a  separate  file  at  the 
beginning  of  the  session. 

2.  Post-processor  on  Disk  2 

The  data  was  processed  using  the  Trimvec  Trim640  program,  Revision 
AB.  Trim640  is  a  relative  positioning,  post-processing  program  that  provides 
triple  and  double  difference  solutions  for  two  sets  of  the  4000SX  carrier  phase 
data  logged  simultaneously  at  two  stations.  Trim640  adopts  the  best  C/A  code 
position  during  the  data  loading.  Only  the  broadcast  ephemeris  can  be  used  to 
compute  fixed  orbit  satellite  positions.  Trim640  limits  processing  to  700 
epoches,  so  the  first  700  epoches  for  each  observing  session  are  used.  The 
antenna  on  the  roof  of  Building  224  on  the  NPS  campus  was  used  as  reference 
station  and  its  coordinates  were  kept  fixed.  Permanent  bench  marks  were  chosen 
as  Trim640  reference  stations  for  the  off  campus  area  because  the  observing 
sessions  at  them  were  shorter  than  the  observing  session  at  Building  224  on  the 
NPS  campus. 

3.  Geoid.exe  on  Disk  3 

This  program  computes  global  geoid  height,  N0(X,Y,Z),  at  any  point 

with  an  accuracy  of  a  few  meters  [Trimble,  1987b].  The  global  geoid  heights, 
N0(X,Y,Z),  are  based  on  Rapp's  1978  360  x  360  model,  an  harmonic  expansion 

of  the  geopotential  referred  to  GRS  80.  Heights  are  suitable  for  showing  relative 
shape  and  trends  in  geoid.  Global  geoid  height  for  the  Monterey  Bay  area  from 
36°  34*  N  to  36°  49'  N  latitude,  and  from  121°  34'  W  to  121°  55'  W  longitude  is 
given  in  Figure  6.  A  Fortran  program  GPSCON  (Appendix  A)  was  used  for  the 
contoured  plot. 

A  second  program  also  named  Geoid.exe  was  obtained  from  the  NGS  of 
NOAA's  National  Ocean  Service  Charting  and  Geodetic  Services  uses  the  same 
Rapp's  1986  360  x  360  model.  The  global  geoid  height  in  the  Monterey  Bay  area 
from  the  NGS  program  is  shown  in  Figure  7.  Because  the  NGS  program  rounds 
off  to  the  nearest  decimeters,  there  are  some  steps  on  the  curves.  The  differences 
between  Figure  6  and  7  may  be  attributed  to  differences  between  Rapp's  1978 
and  1986  models. 
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Figure  6.  The  global  geoid  height  in  the  Monterey  Bay  area  calculated 
using  the  Trimvec  Geoid.exe  program. 
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Figure  7.  The  global  geoid  height  in  the  Monterey  Bay  area  calculated 
using  the  NGS  Geoid.exe  program. 
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4.   Satellite  Visibility  Program  on  Disk  3 

The  program  Stvis.exe  is  a  system  for  generating  visibility  data  for  the 
GPS  satellites.  The  program  takes  as  input  the  observer's  latitude,  longitude,  a 
mnemonic  name  for  the  location,  data  for  calculations  and  optionally  the  time  zone 
offset  from  Universal  Time  Coordinates  (UTC).  NPS  in  Monterey,  California,  was 
used  as  the  reference  station  (  36°  35'  56.1"  N,  121°  52'  36.0"  W,  and  altitude  -19 
m).  Satellite  orbit  data  is  contained  in  the  file,  Almanac.dat,  which  is  produced  by 
processing  data  collected  from  the  actual  satellites.  The  satellite  visibility  chart 
shows  the  tracks  of  the  GPS  satellites  during  the  planned  observation  session 
(Figure  8). 


Figure  8.   Sky  plots  of  satellite  tracks  for  the  Monterey  Bay  area. 

Elevation  angles  are  dotted  concentric  circles.  Zenith  is  at  the  center. 
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V.   DATA   PROCESSING 


A.    ORTHOMETRIC  HEIGHT  ON  NPS  CAMPUS 

In  calculating  the  most  probable  elevations  for  the  temporary  bench  marks, 
station  TREE  was  used  as  a  reference,  where  the  elevation  was  assumed  to  be  0 
meters  above  MSL.  The  interlocking  levelling  circuit  on  the  NPS  campus  is 
shown  in  Figure  9. 


GHIq 

—         H$    TREE 

\GH6 

GH2  <i 

GH7              \ 

GH8 

®  GH5 

GH3 

GH4 

Figure  9.     The  interlocking  levelling  circuit  on  the  NPS  campus. 

A  constrained  least  squares  adjustment  was  used  to  do  the  computation.  A 
weight  of  100  was  assumed  at  the  TREE  and  weights  of  1  were  used  for  all 
observations.  The  13  observation  equations  are: 
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The  observation  equations  in  matrix  form  are  [Wolf,  1980]: 


P  A  H  =    P  L    +    PV 


The  weight  matrix  P13xl3,  coefficient  matrix  A13x9,  observation  matrix  L13xl, 
and  residual  matrix  V13xl  are: 
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The  program  Lobs. basic  [Jeyapalan,  1988]  was  used  to  find  H  for  the  temporary 
bench  marks.  The  program  (Appendix  B)  gives  the  H  and  V  matrices: 
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H  = 


^0.000\ 
-0.808 
4.579 
7.728 
9.262 
8.246 
4.032 
4.521 

V  8.445  7 


V  = 


f  0.00000  ^ 
0.00007 
0.00007 
-0.00022 
-0.00044 
-0.00022 
-0.00022 
0.00007 
-0.00019 
-0.00010 
-0.00010 
0.00022 

V  0.000237 


The  standard  deviation  of  the  observation  equations  is  0.36  mm. 

B.    GPS  DATA  PROCESSING 

1.  Automatic  Processing  Mode 

Trimble  recommends  that  an  automatic  processing  mode  be  used  if  more 
than  one  hour  of  data  from  four  or  five  GPS  satellites  is  taken  at  each  site  for 
baselines  lengths  up  to  30  km. 

The  output  files  contain  triple  difference,  double  difference  float  and 
double  difference  fixed  solutions.  The  double  difference  fixed  solution  provides 
the  most  accurate  solution  if  the  following  conditions  are  met: 

a.  Integer  bias  search  indicates  the  ratio  sum-of- squares  must  be  greater  than 
3.0; 

b.  RMS  (cycles)  less  than  0.05;  and 

c.  Difference  between  the  float  and  fixed  solutions  is  less  than  10  cm,  in  any 
component  of  the  baseline  (X,Y,or  Z). 

In  the  GPS  data  processing  all  the  solutions  from  the  bench  marks 
satisfied  these  conditions  except  station  J  697  for  which  RMS  was  0.068.  The 
distance  from  NPS  Building  224  to  J  697  is  33  km.  Trimble  recommends  a 
baseline  length  greater  than  30  km  when  the  fixed  solution  does  not  meet  the 
above  conditions.  The  float  solution  should  be  used  provided  the  RMS  of  fit  is 
better  than  0.08.  Ah  of  double  difference  fixed  solutions  were  used  for  the  bench 
marks  in  this  study.  The  double  difference  fixed  solutions  are  shown  in  Table  3. 
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Table  3 

.     DOUBLE  DIFFERE 

NCE  FIXED  SOLUTIONS 

Bench 

Slope  distance 

Coordinates  difference  between 

mark 

from  Bldg.  224 

RMS 

Ratio 

fixed  and  float  solution 

(m) 

(cycles) 

AX  (m)      AY  (m)       AZ  (m) 

TREE 

66 

0.039 

20.6 

-2.0 

0.3 

-1.4 

GH1 

98 

0.014 

257.3 

-0.5 

0.4 

0.0 

GH2 

256 

0.025 

94.7 

-3.2 

3.7 

0.3 

GH3 

536 

0.014 

276.4 

-1.1 

0.6 

-0.1 

GH4 

547 

0.019 

151.0 

-0.7 

0.4 

-0.2 

GH5 

424 

0.014 

245.7 

0.0 

-1.6 

-1.2 

GH6 

166 

0.015 

254.9 

-0.7 

0.5 

-0.1 

GH7 

169 

0.015 

288.6 

-0.7 

0.3 

-0.1 

GH8 

413 

0.016 

160.0 

-0.4 

0.5 

0.1 

K152 

14914 

0.035 

9.7 

-0.3 

0.1 

-0.4 

B21 

2082 

0.014 

193.5 

0.4 

-0.2 

0.1 

S812 

17062 

0.020 

15.8 

-1.2 

1.7 

0.1 

P812 

20011 

0.020 

31.7 

0.9 

-1.4 

0.4 

D697 

25249 

0.045 

11.1 

8.2 

9.3 

0.4 

J  697 

32692 

0.068 

6.3 

13.8 

-10.7 

2.6 

F813 

16861 

0.022 

35.9 

7.1 

-5.2 

0.1 

GWM27 

16687 

0.049 

15.3 

4.4 

-5.7 

0.5 

Default  parameters  are  assumed  for  automatic  processing.  These 
parameters  have  a  minimum  elevation  mask  of  15°  for  double  differences. 
Station  one  coordinates  were  the  best  code  positions  during  data  logging.  Actual 
surface  meteorological  values  were  used  in  the  processing. 

For  sets  of  data  covering  more  than  one  hour,  epoch  increments  of  five 
or  ten  provide  full  precision.  Automatic  processing  begins  with  several  iterations 
of  triple  differences  to  establish  a  starting  value  for  the  double  difference 
solution.  A  cycle  slip  fixer  processes  every  epoch. 

Double  difference  solutions  begin  with  several  iterations  of  the  float 
solution  where  biases,  as  well  as  baseline  components,  are  computed.  It  is  called 
the  float  solution,  since  the  biases  are  allowed  to  float  and  are  not  constrained  to 
be  integers.  Each  iteration  should  indicate  convergence  toward  zero  for 
observed  minus  computed  phases.     The  default  elevation  mask  is   15°  for 
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processing  double  differences.  After  completion  of  the  float  solution, 
correlations  and  biases  are  computed.  The  integer  search  algorithm  begins  to 
determine  the  proper  integer  values  for  the  bias.  After  the  biases  are  set  to  their 
integer  values,  the  processing  is  repeated  to  compute  the  fixed  solution.  Ellipsoid 
height  differences  Ah  for  the  temporary  bench  marks  with  respect  to  NPS 
Building  224,  found  using  double  difference  fixed  solutions  is  given  in  Table  4. 


Table     4.     Ah     AND     WGS     84     COORDINATES 
TEMPORARY  BENCH  MARKS. 


FOR     THE 


Bench  mark 

Ah  (m) 

X(m) 

Y(m) 

Z(m) 

TREE 

-8.3532 

-2707298.015 

-4353450.286 

3781781.432 

GH1 

-9.1814 

-2707385.942 

-4353407.659 

3781790.139 

GH2 

-3.7766 

-2707502.139 

-4353548.202 

3781549.491 

GH3 

-0.6382 

-2707528.391 

-4353750.407 

3781314.047 

GH4 

0.9055 

-2707367.536 

-4353813.483 

3781306.605 

GH5 

-0.1099 

-2707202.581 

-4353799.531 

3781493.061 

GH6 

-4.3397 

-2707252.929 

-4353587.524 

3781666.956 

GH7 

-3.8509 

-2707406.088 

-4353553.901 

3781589.428 

GH8 

0.0843 

-2707329.452 

-4353759.387 

3781428.906 

Ah  for  the  permanent  bench  marks  with  respect  to  NPS  Building  224 
using  double  difference  fixed  solutions  are  given  in  Table  5. 
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Table     5.     Ah     AND     WGS     84 
PERMANENT  BENCH  MARKS 


COORDINATES     FOR     THE 


Bench  mark 

Ah  (m) 

X(m) 

Y(m) 

Z(m) 

K152 

2.2398 

-2696890.306 

-4350931.132 

3792065.393 

B21 

-8.4275 

-2708436.436 

-4351973.362 

3782674.213 

S812 

1.6335 

-2692137.111 

-4360870.528 

3784094.277 

P812 

-1.3569 

-2689245.909 

-4360679.825 

3786312.651 

D697 

15.1728 

-2685153.665 

-4357222.820 

3793175.171 

J  697 

71.3023 

-2679276.491 

-4356508.580 

3798216.777 

F813 

-6.7070 

-2695613.779 

-4350454.843 

3793463.222 

GWM27 

49.1051 

-2692408.721 

-4360427.654 

3784433.846 

2.  Obtaining  Geoid  Height 

Geoid  heights  of  bench  marks  were  obtained  by  using  the  Geoid.exe 
program  from  NGS  or  Trimvec  (Table  6). 


Table  6.    GEOID  HEIGHTS 


Bench  mark 

NGS  (m) 

Trimvec  (m) 

TREE 

-34.7 

-34.30428 

GH1 

-34.7 

-34.30671 

GH2 

-34.7 

-34.31508 

GH3 

-34.7 

-34.32072 

GH4 

-34.7 

-34.31605 

GH5 

-34.7 

-34.30745 

GH6 

-34.7 

-34.30532 

GH7 

-34.7 

-34.31140 

GH8 

-34.7 

-34.31247 

K152 

-34.1 

-33.83303 

B21 

-34.7 

-34.31918 

S812 

-33.9 

-33.86473 

P812 

-33.8 

-33.76630 

D697 

-33.6 

-33.58265 

J  697 

-33.3 

-33.42047 

F813 

-34.0 

-33.77899 

GWM27 

-33.9 

-33.86547 
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3.  Determination  of  Local  Geoid  Model 

The  method  of  collocation  was  used  to  solve  the  local  geoid  model.  The 
procedure  is  described  as  follows  : 

a.   Reference  Station  Selection 

Station  TREE  was  selected  as  a  reference  station  on  the  NPS  campus. 
Coordinate  differences  AX,  AY  and  AZ  in  WGS  84,  for  the  temporary  bench 
mark  with  respect  to  the  TREE  were  computed.  These  are  shown  in  Table  7. 


Table  7.    AX,  AY  AND  AZ  ON  THE  NPS  CAMPUS 


Bench  mark 

AX(m) 

AY(m) 

AZ(m) 

TREE 

0.000 

0.000 

0.000 

GH1 

-87.927 

42.627 

8.707 

GH2 

-204.124 

-97.916 

-231.941 

GH3 

-230.376 

-300.121 

-467.385 

GH4 

-69.521 

-363.197 

-474.827 

GH5 

95.434 

-349.245 

-288.371 

GH6 

45.085 

-137.238 

-114.476 

GH7 

-108.073 

-103.615 

-192.004 

GH8 

-31.437 

-309.101 

-352.526 

Station  K  152  was  selected  as  a  reference  station  in  the  large  off 
campus  area.  The  coordinate  differences,  AX,  AY  and  AZ  between  the 
permanent  bench  marks  and  K  152  were  computed  in  WGS  84  (Table  8). 
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Table  8.    AX,  AY  AND  AZ  IN  THE  OFF  CAMPUS  AREA 


Bench  mark 

AX(m) 

AY(m) 

AZ(m) 

K152 

0.000 

0.000 

0.000 

B21 

-11546.130 

-1042.230 

-9391.180 

S812 

4753.195 

-9939.396 

-7971.116 

P812 

7644.397 

-9748.693 

-5752.742 

D697 

11736.641 

-6291.688 

1109.778 

J  697 

17613.815 

-5577.448 

6151.384 

F813 

1276.527 

476.289 

1397.829 

GWM27 

4481.585 

-9496.522 

-7631.547 

b.  Determine  the  Local  Geoid  Model 

I  start  with  assuming  the  signal,  S,  equals  zero  and  the  weight  matrix 
P,  equals  the  unit  matrix.  Seven  control  marks  in  both  areas  led  to  seven 
observation  equations.  The  coordinate  differences  between  the  reference  marks, 
which  are  TREE  on  the  NPS  campus  and  K  1 52  in  the  off  campus  area,  are  the 
coefficients  of  the  parameters.  The  observation  equation  is  given  by 

AN  =  H  +  NQ(X,Y,Z)  -  Ah 

where 

H  :  the  orthometric  height,  from  levelling 

N0(X,Y,Z)  :  the  global  geoid  height,  obtained  form  NGS  or  Trimvec 

Ah  :  the  ellipsoid  height  difference  with  respect  to  the  reference  station 
The  local  variation  of  the  geoid  height,  AN,  is  also  given  by 

AN  =  h   +  combination  of  coordinate  differences 

h0,  the  ellipsoid  height,  includes  a  constant  correction  to  N0(X,Y,Z)  at  the 
reference  station.  AN  for  temporary  bench  marks,  N0(X,Y,Z)  from  the  NGS 
(designated  as  NAN)  and  N0(X,Y,Z)  from  the  Trimvec  (designated  as  TAN)  are 
given  in  Table  9. 


38 


Table  9.    AN  ON  THE  NPS  CAMPUS 

Bench  mark 

NAN  (m) 

TAN  (m) 

TREE 

-26.3468 

-25.9511 

GH1 

-26.3266 

-25.9333 

GH2 

-26.3444 

-25.9595 

GH3 

-26.3338 

-25.9545 

GH4 

-26.3435 

-25.9595 

GH5 

-26.3441 

-25.9515 

GH6 

-26.3283 

-25.9336 

AN  for  permanent  bench  marks  with  N0(X,Y,Z)  from  the  NGS  (designated  as 
NAN)  and  N0(X,Y,Z)  from  the  Trimvec  (designated  as  TAN)  are  given  in  Table 
10. 

Table  10.    AN  IN  THE  OFF  CAMPUS  AREA 


Bench  mark 

NAN  (m) 

TAN  (m) 

K152 

-19.3748 

-19.1078 

B21 

-20.4855 

-20.1047 

S812 

-19.5555 

-19.5202 

P812 

-19.3091 

-19.2754 

D697 

-18.7158 

-18.6984 

J  697 

-19.5413 

-19.6618 

F813 

-19.3100 

-19.0890 

A  least  squares  adjustment  program  Lobs. basic  was  used  to  determine 
the  geoid  model  which  has  the  smallest  standard  deviation.  This  was  done  by 
using  the  different  combinations  of  the  coordinate  differences.  GPS  data  from 
the  off  campus  area  were  studied  to  determine  the  local  geoid  model.  N0(X,Y,Z) 

from  the  Trimvec  program  were  used  to  compute  AN. 
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Four-parameter  combinations  and  their  standard  deviations  for  the  off 
campus  area  are  shown  in  Table  1 1 . 

Table  11.    FOUR  PARAMETERS 


Parameters 

Standard  deviation  (m) 

hQ  +  aAX  +  bAY  +  cAZ 

0.45 

hQ  +  aAX  +  bAY  +  cAY2 

0.52 

hQ  +  aAX  +  bAY  +  cAX2 

0.34 

hQ  +  aAX  +  bAY  +  cAXAY 

0.47 

hQ  +  aAX  +  bAZ  +  cAY2 

0.52 

hQ  +  aAX  +  bAZ  +  cAX2 

0.34 

hQ  +  aAX  +  bAZ  +  cAXAY 

0.47 

hQ  +  aAX  +  bAY2  +  cAX2 

0.31 

hQ  +  aAX  +  bAX2  +  cAXAY 

0.37 

hQ  +  aAY  +  bAZ  +  cAY2 

0.52 

hQ  +  aAY  +  bAZ  +  cAX2 

0.34 

hQ  +  aAY  +  bAZ  +  cAXAY 

0.46 

hQ  +  aAY  +  bAY2  +  cAX2 

0.43 

hQ  +  aAY+  bAY2  +  cAXAY 

0.58 

hQ  +  aAY+  bAX2  +  cAXAY 

0.28 

ho  +  aAZ  +  bAY2  +  cAX2 

0.37 

hQ  +  aAZ  +  bAX2  +  cAXAY 

0.32 

hQ  +  aAZ  +  bAY2  +  cAXAY 

0.49 

hQ  +  aAY2  +  bAX2  +  cAXAY 

0.20    * 

*  The  smallest  standard  deviation. 
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Five-parameter  combinations  and  their  standard  deviations  for  the  off  campus 
area  are  given  in  Table  12. 

Table  12.    FIVE  PARAMETERS 


Parameters 

Standard  deviation  (m) 

hQ  +  aAX  +  bAY  +  cAZ  +  dAX2 

0.43 

hQ  +  aAX  +  bAY  +  cAZ  +  dAY2 

0.57 

hQ  +  aAX  +  bAY  +  cAZ  +  dAXAY 

0.39 

hQ  +  aAX  +  bAZ  +  cAX2  +  dAY2 

0.57 

hQ  +  aAX  +  bAZ  +  cAY2  +  dAXAY 

0.59 

hQ  +  aAX  +  bAX2  +  cAY2  +  dAXAY 

0.22 

hQ  +  aAX  +  bAZ  +  cAX2  +  dAY2 

0.57 

hQ  +  aAY  +  bAZ  +  cAX2  +  dAXAY 

0.34 

hQ  +  aAY  +  bAZ  +  cAX2  +  dAXAY 

0.34 

hQ  +  aAY  +  bAZ  +  cAX2  +  dAY2 

0.56 

hQ  +  aAY  +  bAX2  +  cAY2  +  dAXAY 

0.10    * 

hQ  +  aAZ  +  bAX2  +  cAY2  +  dAXAY 

0.19 

*  The  smallest  standard  deviation. 


Six-parameter  combinations  and  their  standard  deviations  for  the  off  campus  area 
are  given  in  Table  13. 

Table  13.     SIX  PARAMETERS 


Parameters 

Standard  deviation  (m) 

hQ  +  aAX  +  bAY  +  cAZ  +  dAX2  +  eAY2 

0.26 

hQ  +  aAX  +  bAY  +  cAZ  +  dAX2  +  eAXAY 

0.59 

hQ  +  aAX  +  bAY  +  cAX2  +  dAY2  +  eAXAY 

0.13 

hQ  +  aAX  +  bAZ  +  cAX2  +  dAY2  +  eAXAY 

0.12    * 

hQ  +  aAY  +  bAZ  +  cAX2  +  dAY2  +  eAXAY 

0.14 

*  The  smallest  standard  deviation. 
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For  seven  observation  equations  and  seven  parameters  the  degree  of  freedom 
equals  zero,  so  the  standard  deviation  is  undefined.  Seven-parameter 
combinations  and  their  standard  deviations  for  the  off  campus  area  are  given  in 
Table  14. 

Table  14.     SEVEN  PARAMETERS 


Parameters 

Standard  deviation  (m) 

hQ  +  aAX  +  bAY  +  cAZ  +  dAX2  +  eAY2  +  f AXAY 

undefined 

The  best  standard  error  of  the  four-parameter  geoid  model  is  greater 
than  for  five-parameter  or  six-parameter  geoid  models.  The  seven-parameter 
geoid  model  did  not  converge,  so  five-parameter  and  six-parameter  geoid  models 
were  chosen  for  this  study.  The  five-parameter  geoid  model  is  given  by 

N  =  NQ(X,Y,Z)  +  hQ  +  aAY  +  bAX2  +  cAY2  +  dAXAY 

and  the  six-parameter  geoid  height  is  given  by 

N  =  NQ(X,Y,Z)  +  hQ  +  a'AX  +  b'AZ  +  c'AX2  +  d'AY2  +  e'AXAY 

c.   Evaluating  the  Geoid  Model  Using  Check  Marks 

In  order  to  evaluate  the  accuracy  of  the  geoid  model,  the  geoid  height 
model  can  be  rearranged  to  calculate  the  orthometric  height  of  check  marks 
which  are  GH7  and  GH8  on  the  NPS  campus  and  GWM  27  in  the  off  campus 
area. 


or 


H  =  h0  +  Ah  -  N0(X,Y,Z)  +  aAY  +  bAX2  +  cAY2  +  dAXAY 

H  =  h0  +  Ah  -  N0(X,Y,Z)  +  a'AX  +  b'AZ  +  c'AX2  +  d'AY2  +  e'AXAY 


The  h0,  a,  b,  c  and  d  of  the  five-parameter  geoid  model,  and  h0,  a',  b',  c',  d'  and 

e'  of  the  six-parameter  geoid  model,  can  be  solved  by  least  squares  adjustments. 
The  results  for  h0,  a,  b,  c  and  d  for  the  five-parameter  model  using  the  seven 
control  marks  in  the  off  campus  area  are  given  in  Table  15.  The  results  for  h0, 
a,  b,  c  and  d  are  shown  in  the  NGS  column,  where  the  N0(X,Y,Z)  was  from  the 
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NGS  Geoid.exe  program.  The  results  for  h0,  a,  b,  c,  d  and  standard  deviation,  o, 
are  shown  in  the  Trimvec  column,  where  the  N0(X,Y,Z)  was  from  the  Trimvec 
Geoid.exe  program. 


Table  15.    h0,  a, 


b,  c,  d,  a  FOR  SEVEN  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.26895 

-19.00915 

a 

-2.344959E-04 

-2.736598E-04 

b 

-8.523017E-09 

-8.528133E-09 

c 

-3.5743 14E-08 

-3.905052E-08 

d 

-2.211550E-08 

-1.728313E-08 

a 

0.130518 

0.095723 

The  results  of  h0,  a,  b,  c,  d  and  c  of  the  five-parameter  geoid  model 
using  the  six  control  marks  (excluding  B  21)  in  the  off  campus  area  are  given  in 
Table  16. 


Table  16.    h0,  a,  b,  c,  d,  a  FOR  SIX  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.25415 

-19.01411 

a 

-2.752543E-04 

-2.600850E-04 

b 

-7.735253E-09 

-8.792540E-09 

c 

-3.777950E-08 

-3.836067E-08 

d 

-1.798617E-08 

-1.867193E-08 

o 

0.171919 

0.133516 

The  results  of  h0,  a,  b,  c,  d  and  a  of  the  five-parameter  geoid  model 
using  the  six  control  marks  (excluding  J  697)  in  the  off  campus  area  are  given  in 
Table  17. 


43 


Table  17.    h0,  a, 

b,  c,  d,  a  FOR   SIX  CONTROL  MARKS 

Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.265355 

-19.03487 

a 

-2.503395E-04 

-1.977533E-04 

b 

-8.731945E-09 

-7.526410E-09 

c 

-3.711466E-08 

-3.246532E-08 

d 

-2.179058E-08 

-1.882791E-08 

a 

0.1841317 

0.1206984 

The  results  of  h0,  a,  b,  c  and  d  of  the  five-parameter  geoid  model  using 
the  five  control  marks  (excluding  B  21  and  J  697)  in  the  off  campus  area  are 
given  in  Table  18.  The  standard  deviation  is  undefined  since  degree  of  freedom 
equals  zero. 

Table  18.    h0,  a,  b,  c,  d  FOR  FIVE  CONTROL  MARKS    


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.37497 

-19.10791 

a 

1.032352E-04 

3.397465E-05 

b 

8.119969E-09 

3.521564E-09 

c 

7.377821E-09 

-3.339665E-09 

d 

1.338776E-09 

-3.667083E-09 

The  results  of  h0,  a,  b,  c,  d  and  a  of  the  five-parameter  geoid  model 
using  the  seven  control  marks  on  the  NPS  campus  are  given  in  Table  19. 


Table  19.    h0,  a,  b,  c,  d,  o  FOR   SEVEN   CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-26.33543 

-25.94105 

a 

2.976507E-06 

1.192093E-07 

b 

-5.419133E-08 

-1.763692E-07 

c 

-3.601599E-08 

-8.489588E-08 

d 

6.30971 1E-08 

-3.702007E-08 

a 

0.013773 

0.014612 
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The  results  of  h0,  a',  b\  c\  d',  e'  and  a  of  the  six-parameter  geoid  model 
using  the  seven  control  marks  in  the  off  campus  area  are  shown  in  Table  20. 

Table  20.    hQ,  a',  b',  c',  d',  e',  c   FOR   SEVEN   CONTROL   MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.27388 

-19.01660 

a' 

2.083 182E-04 

1.842380E-04 

b* 

-2.5615 10E-04 

-2.520234E-04 

c' 

-7.731387E-09 

-8.335974E-09 

d' 

-3.767036E-08 

-3.960304E-08 

e' 

-1.249646E-08 

-1.521767E-08 

a 

0.137149 

0.123960 

The  results  of  h0,  a',  b',  c',  d'  and  e'  of  the  six-parameter  geoid  model 
using  the  six  control  marks  (excluding  B  21)  in  the  off  campus  area  are  given  in 
Table  21.  The  standard  deviation  is  undefined  since  degree  of  freedom  equals 
zero. 

Table  21.    h0,  a',  b',  c',  d',  e'  FOR  SIX  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.37489 

-19.10785 

a' 

2.57641 1E-04 

2.288222E-04 

b' 

-1.691282E-04 

-1.734756E-04 

c' 

-1.083072E-08 

-1.113654E-08 

d' 

-2.817797E-08 

-3.102741E-08 

e' 

-5.824404E-09 

-9.185896E-09 

The  results  of  h0,  a',  b',  c',  d'  and  e'  of  the  six-parameter  geoid  model 
using  the  six  control  marks  (excluding  J  697)  in  the  off  campus  area  are  given  in 
Table  22.  The  standard  deviation  is  undefined  since  degree  of  freedom  equals 
zero. 
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Table  22.    h0,  a',  b',  c',  d',  e'  FOR  SIX  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-19.37466 

-19.10774 

a' 

1.065433E-04 

9.238720E-05 

b' 

-4.604459E-05 

-6.231666E-05 

c' 

-2.012712E-09 

-3.170499E-09 

d' 

-1.141234E-08 

-1.586523E-08 

e' 

-2.528395E-09 

-6.206392E-09 

The  results  of  h0,  a',  b',  c',  d',  e'  and  a  of  the  six-parameter  geoid  model 
using  the  seven  control  marks  on  the  NPS  campus  are  given  in  Table  23. 

Table  23.    hft,  a',  b',  c',  d\  e',  a   FOR   SEVEN  CONTROL  MARKS 


Parameter 

NGS  (m) 

Trimvec  (m) 

ho 

-26.33389 

-25.93820 

a' 

1.645088E-05 

4.556775E-05 

b' 

4.225  970E-05 

6.1 6908 1E-05 

c' 

6.856863E-08 

6.391201E-08 

d* 

6.00703 1E-08 

5.878974E-08 

e' 

1.781 154E-07 

1.767185E-07 

a 

0.018859 

0.018872 

Fortran  program  GPSDIS  (Appendix  C)  was  used  to  calculate  the 
orthometric  heights  of  the  check  marks.  The  results  are  discussed  in  the  Chapter 
VII. 

d.   Obtaining  the  Weight  Matrix 

Since  the  standard  error  of  the  geoid  model  is  about  ±  2  cm,  which  is 
the  same  or  better  than  GPS  Ah  accuracy,  it  is  not  actually  necessary  to  do 
further  computation.  Nevertheless,  the  weight  matrices,  P,  were  computed. 

The  five-parameter  geoid  model  used  for  the  NPS  campus  and  the 
six-parameter  geoid  model  used  for  the  off  campus  area  were  employed  in 
obtaining  P.  The  covariance  function,  C(r),  of  the  control  marks  (p.  18)  can  be 
obtained  by  using  the  least  squares  residuals.  The  coefficients  B  and  k  of  C(r) 
can  be  obtained  by  solving  simultaneously  the  equations  for  C(r)  with  different 
distances  from  the  centers. 
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(1)  P  for  NPS  Campus.  The  residuals  of  the  five-parameter  geoid 
model  of  the  control  marks  on  the  NPS  campus  are  given  in  Table  24.  The 
global  geoid  heights  were  from  the  NGS. 

Table  24.     RESIDUALS  FOR  FIVE-PARAMETER  MODEL 


Control  mark 

#  of  residual 

Residual  (m) 

TREE 

v, 

1.137352E-02 

GH1 

v2 

-9.420395E-03 

GH2 

v3 

7.339478E-03 

GH3 

V4 

-4.278183E-03 

GH4 

v5 

3.572464E-03 

GH5 

v6 

6.446839E-04 

GH6 

V7 

-8.714676E-03 

For  r   =  593  m 


ViV,  +  vv  +  wwA  +  vyK  +  v  v 

CCfj)  =    — L- ^ — *r-* ^^  =-2.194272E-03 


For  r   =  319  m 

2 


v,v +V  V+V  V+V  V +V  V +V  V +V  V+V  V +V  V, 
C(r2)=     ]    3      2   3      2   7      3   4      3g  5      3    6      3   7     4   6      5    6  =  3.945813E-6 


Thus  the  equations  of  covariance  can  be  written  as 


and 


C(rj)  =  0.0138  +  B  sinCkr^ 


C(r2)  =  0.0138  +  B  sin(kr2) 


Approximate  solutions  for  coefficients  B  and  k  can  be  found  by  using  the 
Maclaurin  series  expansion.  The  solved  covariance  function  is  given  by 

C(r)  =0.0138- 0.0160  sin(1.7670r) 
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C    can  be  obtained  by  using  the  C(r)  with  the  distances  between  the  control 

q  ° 

marks.  Fortran  program  DISTCO  (Appendix  D)  was  used  to  compute  the 
distances  between  the  control  marks,  as  well  as  C  and,  hence,  the  weight  matrix, 
P  =  C  _1.  The  Matlab  program  on  the  NPS  mainframe  computer  was  used  to  find 
the  inverse  of  C  (P  =  C  _1).  P  for  the  five-parameter  geoid  model  for  the  NPS 
campus  is  given  by 


s 


P    = 


\ 


162 


ZEROS 


142 


78 


98 


102 


ZEROS 


87 


100 


v 


(2)  P  for  off  Campus  Area.  The  residuals  of  the  six-parameter  geoid 
model  of  the  control  marks  in  the  off  campus  area  are  given  in  Table  25.  The 
global  geoid  heights  were  calculated  using  the  Trimvec  program. 

Table  25.     RESIDUALS  FOR  SIX-PARAMETER  MODEL 


Control  mark 

#  of  residual 

Residual  (m) 

K152 

v, 

9.120178E-02 

B21 

v2 

-9.773254E-03 

S812 

V3 

6.391526E-03 

P812 

v4 

1.983643E-04 

D697 

v5 

-2.779961E-02 

J  697 

v6 

1.685524E-02 

F813 

v, 

-7.65171 1E-02 

For  r  =  13994  m 


vtvo  +  V  V    +  V  V    +  V  V    +  V  v   +  v  v 
C^)  =  -1-2 ^ l    4  5     l    5 3— — 7-  =  -6.659883E-04 
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For  r   =  19230  m 

2 


V,V +V  V+V  V 
C(r2)  =  -L-S — |-J — —    =  7.912463E-04 


Thus  the  equations  of  covariance  can  be  written  as 


and 


C(r  )  =  0.0124  +  B  sinCkr^ 


C(r2)  =  0.0124 +  Bsin(kr2) 


Approximate  solutions  for  coefficients  B  and  k  can  be  found  by  using  a 
Maclaurin  series  expansion.  The  solved  covariance  function  is  given  by 

C(r)   =0.0124-  0.1329  sin(2.8278  r) 

Similarly,  C   can  be  obtained  by  using  the  C(r)  with  the  distances  between  the 

control  marks.  Fortran  program  DISTCO  (Appendix  D)  was  used  to  compute 
the  distances  between  the  control  marks  and  C  ,  to  obtain  the  weight  matrix  (P  = 
C  _1)-  The  Matlab  program  on  the  NPS  mainframe  computer  was  used  to  find  the 
inverse  of  C  (P  =  C  _1).  P  for  the  six-parameter  geoid  model  of  the  off  campus 
area  is  given  by 


s 


P     = 


\ 


25 


^ 


15 


ZEROS 


15 


16 


10 
ZEROS 


25 
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VI.  EVALUATION  OF  GPS  RECEIVER  ERRORS 


Factors  affecting  the  GPS  measured  ellipsoid  height  differences,  Ah, 
include:  (1)  comparing  GPS  differencing  solutions,  (2)  standard  error  of  GPS 
observations,  (3)  corrections  for  surface  meteorological  values,  and  (4) 
observation  durations  for  GPS. 

A.    COMPARING  GPS  DIFFERENCING  SOLUTIONS 

1.  Ellipsoid  Height  Differences  Tested  at  a  Fixed  Position 

In  order  to  determine  the  best  solution  for  Ah  using  GPS,  the  Ah  data 
from  a  Trimble  4000SX  GPS  receiver  were  collected  in  the  K  parking  lot  on  the 
NPS  campus  from  0633  to  1007  UTC  on  February  19,  1988.  The  design  of  the 
experiments  is  shown  in  the  Figure  10. 


Antenna 
—A 


[    Grid  PC 

Trimble      4000SX 

i                              i 

A     A     A    A 

y-i  ¥-»<-!        ¥ 

ZZ  GPS  Locator 

Plumb  bob 


-/- 


1 2-foot  ladder 


Board 


Fixed  point  on  the  ground 


Figure    10.   Illustration   of  the   Ah   test   of  Trimble   4000SX   GPS 
receiver. 
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A  plumb  bob  was  used  to  point  to  a  fixed  position  (36°  35'  56"  N,  121° 
52'  36"  W)  on  the  ground.  The  observations  were  taken  by  changing  the  antenna 
height  over  a  distance  of  about  60  cm,  which  was  measured  with  a  tape  every 
half  hour.  The  ellipsoid  height  differences,  Ah,  and  the  orthometric  height 
differences,  AH,  should  agree  with  the  change  of  antenna  height  for  a  fixed 
position.  Triple  difference  solutions,  AhTRI,  are  given  in  Table  26.  Table  also 
gives  the  differences  between  the  AhTRI  at  succeeding  observation  times,  DAhTRI; 
double  difference  float  solutions,  AhFLT;  the  differences  between  the  AhFLT  at 
succeeding  observation  times,  DAhFLT;  and  double  difference  fixed  solutions, 
AhFK,  and  the  differences  between  the  AhFIX  at  succeeding  observation  times, 

DAhFDC 


Table  26. 

Ah  AND 

DAh  FOR  DIFFERENT  ANTENNA  HEIGHTS 

TIME  (UTC) 

AhTRI 

DAhTRI 

AhpLT 

DAhFLT 

Ah^ 

DAhFK 

(19-2-1988) 

(m) 

(m) 

(m) 

(m) 

(m) 

(m) 

0633  -  0712 

4.4410 

-0.7189 

4.4429 

-0.6339 

4.5045 

-0.5463 

0736  -  0808 

5.1599 

-0.3973 

5.0768 

-0.5676 

5.0508 

-0.5964 

0824  -  0855 

5.5572 

-0.1954 

5.6444 

-0.7144 

5.6472 

-0.7030 

0911  -0943 

5.7526 

-0.4531 

6.3588 

-0.2920 

6.3502 

-0.2660 

0945  -  1007 

6.2057 

6.6508 

6.6162 

The  differences  between  the  GPS  DAh  and  the  differences  of  the 
measured  ellipsoid  height  differences,  DAhMEA,  are  given  in  Table  27. 
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Table  27. 

THE  DIFFERENCE  OF  Ah 

ANTENNA 

DAhMEA 

DAhMEA-DAhTRI 

DAhMEA  -  DAhFLT 

DAhMEA  "  DAhFK 

HEIGHT  (m) 

(m) 

(m) 

Cm) 

(m) 

2.809 

-0.5470 

0.1719 

0.0869 

-0.0007 

2.262 

-0.5960 

-0.1987 

-0.0284 

0.0004 

1.666 

-0.5920 

-0.3966 

0.1224 

0.1110 

1.074 

-0.5900 

-0.1369 

-0.2980 

-0.3240 

0.484 

The  results  show  that  the  double  difference  fixed  solution  is  the  best. 
The  difference  between  the  measured  ellipsoid  height  difference  and  the  GPS 
ellipsoid  height  difference  was  about  1  mm.  The  last  two  results  were  not  good 
because  they  had  multipath  effects  and  imaging  effects  from  the  ground  (antenna 
heights  1.074  m  and  0.484  m),  so  the  solutions  of  AH  compared  to  the  true 
differences  were  large.  This  indicates  that  the  antenna  should  be  at  least  1  m 
above  the  ground. 

2.  Ellipsoid  Height  Differences  Tested  at  Different  Positions 

Since  the  geoid  on  the  NPS  campus  has  a  flat  geoid  slope,  the  differences 
between  Ah,  DAh,  and  the  orthometric  height  difference,  AH,  of  the  temporary 
bench  marks  can  be  neglected.  The  Ah  and  DAh  from  the  GPS  observations  are 
given  in  Table  28.  The  differences  between  the  GPS  DAh  and  the  AH  are  given 
in  Table  29. 
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Table  28.  Ah  AND  DAh  ON  THE  NPS  CAMPUS 

Bench 

AhTRl 

DAhTRI 

Ahpur 

DAhFLT 

Ahpjx 

DAhFIX 

mark 

(m) 

(m) 

(m) 

(m) 

(m) 

(m) 

TREE 

-8.3358 

0.8501 

-8.3375 

0.8313 

-8.3532 

0.8282 

GH1 

-9.1859 

-5.3472 

-9.1688 

-5.3690 

-9.1814 

-5.4048 

GH2 

-3.8387 

-3.1994 

-3.7998 

-3.1624 

-3.7766 

-3.1384 

GH3 

-0.6393 

-1.5241 

-0.6374 

-1.5441 

-0.6382 

-1.5437 

GH4 

0.8848 

0.9777 

0.9067 

1.0171 

0.9055 

1.0154 

GH5 

-0.0929 

4.2655 

-0.1104 

4.2303 

-0.1099 

4.2298 

GH6 

-4.3584 

-0.4932 

-4.3407 

-0.4916 

-4.3397 

-0.4888 

GH7 

-3.8652 

-3.9412 

-3.8491 

-3.9360 

-3.8509 

-3.9352 

GH8 

0.0760 

0.0869 

0.0843 

Table  29.     DIFFERENCES 

BETWEEN  DAh  AND  AH 

Bench 

H 

AH 

AH-DAh-rRj 

AH  -  DAn^ 

AH  -  DAhj^ 

mark 

(m) 

(m) 

(m) 

(m) 

(m) 

TREE 

0.000 

0.808 

-0.0423 

-0.0235 

-0.0204 

GH1 

-0.808 

-5.387 

-0.0398 

-0.0180 

0.0178 

GH2 

4.579 

-3.149 

0.0504 

0.0134 

-0.0106 

GH3 

7.728 

-1.534 

-0.0099 

0.0101 

0.0097 

GH4 

9.262 

1.016 

0.0383 

-0.0011 

0.0006 

GH5 

8.246 

4.214 

-0.0515 

-0.0163 

-0.0158 

GH6 

4.032 

-0.489 

0.0042 

0.0026 

-0.0002 

GH7 

4.521 

-3.924 

0.0172 

0.0120 

0.0112 

GH8 

8.445 
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It  is  clear  that  the  double  difference  fixed  solutions  give  the  best  solution 
in  a  flat  geoid  slope  area. 

B.  STANDARD  ERROR  OF  GPS  OBSERVATIONS 

To  determine  the  accuracy  of  GPS,  measurements  were  taken  at  bench  mark 
TREE  on  the  NPS  campus.  Three  measurements,  each  of  about  three  hours 
duration,  were  taken  for  this  analysis.  Double  difference  fixed  solutions  were 
used  to  analyze  the  data.  The  double  difference  fixed  solutions  for  these 
measurements  are  given  in  Table  30. 

Table  30.  Ah        AT  STATION  TREE 

r  IX 


Date 

Time  (UTC) 

AhFDC  (m) 

Feb.    5,  88 
Feb.    7,  88 
Feb.  11,  88 

0730-  1025 
0721  -  1004 
0705  -  0959 

-8.4262 
-8.3532 
-8.3602 

The  standard  error  of  the  mean  of  Ah  (a  /  Vn)  is  about  ±  2  cm  at  TREE  and 
the  standard  error  (a)  of  the  single  measurement  is  about  ±  4  cm  for  three  hours 
of  observation. 

Results  obtained  by  Strange  in  1985  for  a  project  southeast  of  Phoenix, 
Arizona,  show  that  using  multiple  GPS  observations  at  the  same  point  lead  to 
uncertainties  typically  less  than  2  cm  in  Z-direction  over  20-km  distances 
[Zilkoski,  1988]. 

Thus  the  accuracy  of  the  ellipsoid  differences  obtainable  by  GPS 
measurements  is  expected  to  be  about  ±  2  cm  for  about  three  hours  of 
observation. 


C.    CORRECTIONS     FOR     SURFACE    METEOROLOGICAL 

VALUES 

The  meteorological  correction  is  a  function  of  the  atmospheric  refractivity 
as  computed  by  the  surface  meteorological  values  of  pressure,  temperature, 
humidity  and  the  elevation  angle  of  the  satellite.  Larger  corrections  are  required 
for  low  elevation  angles,  as  the  signal  travels  a  longer  path  through  the 
troposphere.    A  modified  Hopfield  troposphere  model  [Fell,  1975]  is  used  for 
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automatic  processing  of  Trim640  software.  The  model  corrects  for  at  least  90% 
of  the  tropospheric  delay  [Remondi,  1984].  Default  surface  meteorological 
values  for  automatic  processing  are  1010  millibars,  20°  C  and  50  %  relative 
humidity.  These  parameters  can  be  reset  before  the  processing.  Trim640  allows 
only  one  pressure,  temperature  and  humidity  entry  for  each  site  per  session. 
Ellipsoid  height  differences,  Ahnx,  were  calculated  for  four  marks  on  the  NPS 

campus  and  two  off  campus  using  the  default  meteorological  parameters  and 
using  the  true  values.  There  are  given  in  Table  31  along  with  differences  found 
between  the  two  methods. 

Table    31.    COMPARISON    OF    SURFACE    METEOROLOGICAL 
CORRECTION 


Bench  mark 

Default  Ahp^  (m) 

True  AhFK  (m) 

Difference  (mm) 

GH2 

-3.7777 

-3.7766 

1.1 

GH3 

-0.6376 

-0.6383 

0.6 

GH4 

0.9070 

0.9055 

1.5 

GH7 

-3.8524 

-3.8509 

1.5 

K152 

2.2398 

2.2430 

-3.2 

J  697 

-71.3179 

-71.3023 

15.6 

The  difference  is  about  ±  2  mm  on  the  NPS  campus  and  about  ±  2  cm  in  the 
Monterey  Bay  area.  The  Trim640  program  contains  no  ionospheric  model.  For 
first-order  relative  positioning  (1  part  in  105)  and  a  baseline  shorter  than  50  km 
ionospheric  delay  can  be  neglected  [Remondi,  1984].  In  this  study  surface 
meteorological  corrections  were  made  during  data  processing  to  get  the  best  Ah 
solutions. 

D.    OBSERVATION  DURATIONS  FOR  GPS 

To  find  out  how  long  observations  must  be  made  to  obtain  the  best  solutions 
in  the  field  when  batteries  are  used  for  power,  the  GPS  data  from  S  812,  J  697 
and  K  152  were  segmented  into  averaging  periods  of  length  nAt,  where  At  =  10 
minutes  and  n  =  1,  2,  3,  ...,10.  Default  meteorological  values  were  used  for  this 
processing.    AhFK  for  various  observation  durations  for  the  stations  is  given  in 

Table  32  and  plotted  in  Figure  1 1  to  Figure  13. 
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Table  32.  Ah,  AS  A  FUNCTION  OF  OBSERVATION  DURATION 


FIX 


Time  (min) 

S812(m) 

J  697  (m) 

K  152  (m) 

10 

1.5309 

71.5291 

2.3135 

20 

1.6141 

70.9697 

2.3102 

30 

1.6441 

71.2185 

2.3062 

40 

1.6539 

71.4514 

2.3065 

50 

1.6558 

71.4762 

2.3060 

60 

1.6498 

71.5343 

2.2563 

70 

1.6464 

71.5309 

2.2504 

80 

1.6349 

71.3407 

2.2504 

90 

71.3179 

2.2548 

100 

2.2430 

1.70 

1.65 

- 
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Figure  11.  Ah  vs.  observation  durations  at  station  S  812. 
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Figure  12.  Ah  vs.  observation  durations  at  station  J  697. 
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Figure  13.  Ah  vs.  observation  durations  at  station  K  152. 


It  is  clear  that  Ah  was  affected  by  the  observation  duration.  Ah  varied  quite 
a  lot  in  the  first  forty  minutes.  After  about  forty  minutes,  Ah  tended  to  close  to 
the  mean  value  of  the  observations.  After  sixty  minutes,  Ah  didn't  vary  much. 
Trimble  recommends,  to  ensure  sufficient  quality  data  and  to  obtain  first-order 
results  on  single  observations,  at  least  one  hour  should  be  taken  for  a  four  or  five 
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satellite  session.  For  baseline  lengths  up  to  30  km  the  conditions  for  the  fixed 
double  solution  of  automatic  processing  will  be  met  if  data  are  collected  for 
ninety  minutes  on  four  or  five  satellites  [Trimble,  1987b].  The  results  shows  that 
the  best  solutions  can  be  improved  with  longer  observation  times.  Normally,  at 
least  sixty  minutes  of  observations  are  required  for  vertical  control. 
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VII.    RESULTS    AND   DISCUSSION 


A.    EVALUATATION  OF  GEOID  MODEL 

To  find  the  accuracy  of  the  geoid  model  the  orthometric  heights  of  the 
check  marks  from  GPS,  H^g,  were  calculated,  and  these  were  compared  to  the 
orthometric  height  of  the  check  marks  from  the  levelling,  HLEVELLING 

I  started  with  the  assumptions  that  C  =  I  and  C  =  0.  Five-parameter  and 
six-parameter  geoid  models  were  used  to  calculate  the  orthometric  height  of  the 
check  marks.  G WM  27  is  a  check  mark  in  the  off  campus  area.  GH7  and  GH8 
are  the  check  marks  on  the  NPS  campus.   N0(X,Y,Z)  was  calculated  using  both 

the  NGS  and  the  Trimvec  programs. 

1.  I^pc-  from  Five-Parameter  Geoid  Model 
a.   Hgpg  in  the  off  Campus  Area 

HGPS  in  GWM  27  was  calculated  by  using  seven  control  marks,  six 
control  marks  (excluding  B  21  or  J  697)  and  five  control  marks  (excluding  B  21 
and  J  697).  Comparisons  of  H  in  GWM  27  are  given  in  Table  33. 

Table    33.    COMPARISONS    OF    H    USING    FIVE-PARAMETER 
MODEL  AT    GWM  27 


GWM  27 

"GPS  "  "LEVELLING     (Cm) 

#  of  control  marks 

Trimvec 

NGS 

7 

-11.55 

-2.24 

6  (excluding  B21) 
6  (excluding  J  697) 

-9.77 
-4.87 

-8.07 
-3.65 

5 

2.14 

6.58 

b.   Hqps  on  the  NPS  Campus 

Hqps  in  GH7  and  GH8  were  calculated  by  using  seven  control  marks. 
Comparisons  of  H  in  GH7  and  GH8  are  given  in  Table  34. 
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Table    34.    COMPARISONS    OF    H    USING    FIVE-PARAMETER 
MODEL  AT  GH7  AND  GH8 


Five  parameters 

"GPS  "  "LEVELLING     (Cm) 

Check  mark 

Trimvec 

NGS 

GH7 

GH8 

0.52 
-0.20 

0.82 
0.00 

2.  Hqps  from  Six-Parameter  Geoid  Model 
a.   Hqps  in  the  off  Campus  Area 

HGPS  in  GWM  27  were  calculated  by  using  seven  control  marks,  and 
six  control  marks  (excluding  B  21  or  J  697).  Comparisons  of  H  at  GWM  27  are 
given  in  Table  35. 

Table    35.    COMPARISONS    OF    H    USING    SIX-PARAMETER 
MODEL  AT    GWM  27 


GWM  27 

"GPS  "  "LEVELLING     (Cm) 

#  of  control  marks 

Trimvec 

NGS 

7 

-12.40 

-11.13 

6  (excluding  B21) 
6  (excluding  J  697 

-9.37 
-3.49 

-7.72 
-1.01 

b.  Hqps  on  the  NPS  Campus 

Hqps  of  GH7  and  GH8  were  calculated  by  using  seven  control  marks. 
Comparisons  of  H  at  GH7  and  GH8  are  given  in  Table  36. 


Table    36.    COMPARISONS    OF 
MODEL  AT  GH7  AND  GH8 


H    USING    SIX-PARAMETER 


Six  parameters 

^GPS  "  ^LEVELLING     (Cm) 

Check  mark 

Trimvec 

NGS 

GH7 
GH8 

1.22 
0.22 

1.22 
0.25 

Since  the  accuracy  of  predicted  geoid  height  at  check  marks  (±  0.2  to  2 
cm)  is  less  than  the  accuracy  of  GPS  observations  (±  2  to  4  cm),  it  may  be 
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concluded  that  P  =  I  is  a  satisfactory  covariance  matrix  and  further  computation 
is  unnecessary. 

B.  ACCURACY 

The  best  accuracy  of  the  five-parameter  geoid  model  is  ±  2  cm  in  the  off 
campus  area  and  is  ±  2  to  10  mm  for  the  NPS  campus.  The  best  accuracy  of  six- 
parameter  geoid  model  is  ±  1  cm  in  the  off  campus  area  and  is  ±  2  to  10  mm  on 
the  NPS  campus. 

C.  DISCUSSION 

The  errors  in  the  geoid  model  include  errors  associated  with  GPS  and  with 
levelling. 

1.  Errors  Associated  With  GPS 

a.  Atmospheric  Delays 

The  time  delay  of  RF  signals  passing  through  the  ionosphere  is  due  to 
a  reduction  in  speed  and  the  bending  of  the  ray,  both  effects  being  due  to 
refraction.  The  overall  delay  in  the  signal  is  nearly  inversely  proportional  to  the 
square  of  the  frequency.  The  ionospheric  delay  may  be  calculated  by  using  two- 
frequency  receivers  (C/A  code,  P  code),  and  corrected  for  with  a  satisfactory 
degree  of  accuracy.  Tropospheric  delays  are  independent  of  frequency.  They 
are  relatively  small  and  can  be  modelled  using  the  elevation  angle  of  the  satellite. 
The  Trimble  4000SX  receiver  only  observes  the  C/A  code  .  The  Trim640 
program  contains  no  ionospheric  modelling.  The  modified  Hopfield  model  [Fell, 
1975]  is  used  for  the  tropospheric  delay.  The  combined  effects  of  unmodelled 
ionospheric  and  tropospheric  errors  are  estimated  to  result  in  a  satellite-to-user 
range  error  of  from  2.4  to  5.2  m  [Milliken,  1980]. 

b.  Selection  of  Appropriate  Control  Marks 

Only  eight  permanent  bench  marks  were  recovered  in  the  off  campus 
study  area.  D  697  was  obstructed  by  a  tree.  F  813  was  set  near  a  steel  bridge. 
Although  offset  levelling  was  performed,  the  ellipsoid  height  could  have  been 
affected  by  offsetting.  The  shape  of  the  bench  marks  were  spread  out  roughly  in 
a  rhomboid  shape.  B  21  and  J  697  are  on  the  ends  of  the  rhomboid.  The 
distance  is  about  33  km  from  B  21  to  J  697.  The  elevation  changes  form  5.787  to 
85.061  m.  If  more  bench  marks  could  be  recovered  and  be  connected  to  these 
marks,  the  accuracy  of  the  seven  control  marks  could  be  improved. 
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c.  Error  of  Antenna  Height  Measurements 

The  slant  heights  of  GPS  antenna  were  measured  by  using  a  steel  tape 
before  and  after  data  collection.    Thus,  the  heights  of  antenna  were  calculated. 
The  systematic  error  of  the  tape,  reading  error,  calculation  error  and  a  slack  tape 
doing  measurement  will  affect  the  accuracy  of  the  ellipsoid  height  difference. 
The  error  of  antenna  height  measurement  is  ±  2  to  10  mm. 

d.  Selecting  observation  period  with  optimal  satellite  geometry 
Because  there  are  only  seven  GPS  satellites  in  operation  at  present, 

the  observing  window  is  only  four  to  five  hours  a  day  for  four  or  five  satellites. 
The  observation  schedule  should  be  planed  early  and  considered  the  travel  time, 
and  time  to  reoccupy  the  marks.  In  this  study  most  of  the  GPS  measurements 
were  made  at  night.  PDOP  were  not  always  optimum. 

e.  Error  in  the  Computing  of  Trim640  Program 

There  is  no  standard  specification  of  the  GPS  software  at  present. 
Round  off  error  in  the  data  processing,  the  data  record  frequency,  and  the 
precision  of  the  ephemeris  could  affect  the  components  of  the  baseline. 
Specifications  of  the  GPS  software  certainly  will  be  developed  in  the  future. 
2.  Errors  in  Levelling 

The  error  in  levelling  would  cause  the  errors  of  the  observation 
equations  when  determining  the  local  geoid  model.  The  maximum  allowable 
closing  error  between  the  forward  and  backward  running  of  a  section  is  3.0  mm 
Vk  for  first-order  class  I  levelling,  and  9.0  mm  Vk  for  third-order  class  I 
levelling,  where  k  is  the  length  of  the  section  measured  in  km  [Bodnar,  1975]. 
The  average  distance  from  K  152  to  the  permanent  bench  marks  is  about  13  km 
in  the  off  campus  study  area.  Thus,  for  first-order  class  I  levelling,  the 
maximum  allowable  closing  error  would  be  about  1 1  mm  in  the  off  campus  area. 
The  first-order  bench  marks  were  set  around  1933.  The  elevations  of  the  marks 
may  have  changed  caused  by  subsidences  or  earthquakes.  Since  the  time  did  not 
permit  rerunning  the  levelling  for  the  off  campus  area,  published  elevations  were 
assumed  correct.  Errors  in  levelling  could  contribute  to  errors  in  the  geoid 
model. 

The  distance  of  the  interlocking  levelling  circuit  on  the  NPS  campus  is 
about  0.8  km.  For  third-order  class  I  levelling,  the  maximum  allowable  closing 
error  is  about  7  mm  on  the  NPS  campus.   The  5  mm  closing  error  was  obtained 
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for  that  region.   The  geoid  height  model  thus  had  a  better  accuracy  on  the  NPS 
campus. 
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VIII.   CONCLUSIONS  AND  RECOMMENDATIONS 


Local  geoid  models  for  the  Monterey,  California,  area  were  determined  by 
GPS  differential  positioning  using  the  method  of  collocation.  The  local  geoid 
models  were  based  on  Rapp's  360  degree  x  360  order  global  geoid  models 
determined  from  gravity  measurements.  Control  data  were  adjusted  by  least 
squares  to  solve  for  the  parameters  in  the  local  geoid  model.  Also  studied  were 
factors  that  affected  the  GPS-measured  ellipsoid  height  differences.  These  included 
(1)  comparing  GPS  differencing  solutions,  (2)  standard  error  of  GPS  observations, 
(3)  corrections  for  surface  meteorological  values,  and  (4)  observation  durations  for 
GPS.  Conclusions  are  the  following: 

1 .  The  accuracy  of  a  local  geoid  height,  N,  is  the  same,  whether  the  NGS  or 
the  Trimvec  version  of  Rapp's  global  geoid  model  is  used,  although  N0(X,Y,Z) 

calculated  using  Trimble's  Trimvec  program  gives  the  best  results  for  points  on  the 
NPS  campus,  and  N0(X,Y,Z)  from  the  NGS  program  gives  the  best  results  for  the 

larger  study  area.  For  practical  use  the  orthometric  heights  can  be  calculated  by 
utilizing  the  local  geoid  model. 

2.  The  areas  where  the  geoid  is  smoothest  lead  to  the  best  results.  This  was 
found  by  comparing  the  local  geoid  model  used  in  the  larger  Monterey  Bay  area  to 
the  local  geoid  model  used  on  the  NPS  campus.  Also  the  density  of  control  marks  on 
the  NPS  campus  was  much  higher  than  for  the  larger  study  area.  If  a  higher 
accuracy  is  required  in  a  large  area  where  the  geoid  may  be  undulating,  more 
control  marks  are  required.  These  control  marks  should  be  strategically  located 
throughout  the  network  in  order  to  determine  the  geoid  slope. 

3.  The  local  geoid  model  does  improve  the  accuracy  of  Rapp's  global  geoid 
model  locally.  The  accuracy  for  Rapp's  model  is  ±  2  m  in  an  area  55  km  x  55  km. 
The  accuracy  for  the  local  geoid  model  is  ±  2  cm  in  an  area  30  km  x  30  km  and  ±  1 
cm  in  an  area  0.5  km  x  0.5  km. 

4.  A  GPS  double  difference  fixed  solution  gives  the  best  ellipsoid  height 
differences,  Ah,  for  baseline  lengths  up  to  30  km. 

5.  The  standard  error  of  GPS  observations  of  the  ellipsoid  height  difference, 
Ah,  is  ±  2  cm  to  ±  4  cm,  so  there  is  no  need  to  perform  the  collocation  with  a  non- 
unit  weight  matrix.  Also  it  is  not  possible  to  reduce  the  error  of  the  GPS 
observation  below  ±  2  cm. 
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6.  There  is  no  need  for  local  surface  meteorological  corrections  to  the  GPS 
observations.  Default  meteorological  values  (1010  millibars,  20°  C  and  50% 
relative  humidity)  are  sufficient  to  meet  the  required  accuracy. 

7.  The  longer  the  duration  of  the  GPS  measurements,  the  better  the  result  will 
be.  The  accuracy  of  the  geoid  model  for  the  NPS  campus,  where  GPS 
measurements  were  made  over  four-  to  five-hour  periods,  is  better  than  the 
accuracy  of  the  geoid  height  model  in  the  larger  study  area,  where  measurements 
were  made  for  periods  of  only  ninety  minutes. 

8.  A  five-parameter  geoid  model  gives  the  best  solution  for  the  NPS  campus, 
while  a  six-parameter  geoid  model  gives  the  best  solution  for  the  larger  study  area. 
The  five-parameter  geoid  model  which  does  not  include  a  AZ  term  should  be  used 
with  caution  where  AZ  greatly  exceeds  5  m. 

Recommendations  are  as  follows: 

1 .  The  GPS  antenna  should  be  set  at  least  1  m  above  the  ground  to  reduce  the 
multiple  effects  from  the  ground.  The  height  of  the  antenna  should  be  carefully 
measured  before  and  after  the  data  collection. 

2.  The  configuration  of  the  control  marks  should  be  carefully  planned.  Sites 
should  be  selected  to  optimize  connections  to  stations  with  known  precise 
orthometric  height  differences,  minimizing  the  effects  of  obstructions.  Control 
marks  in  a  large  area  should  be  numerous  enough  to  meet  the  required  accuracy. 

3.  To  minimize  computation  errors  least  squares  adjustments  should  be  made 
using  double  precision  processing  when  solving  the  geoid  model.  For  large 
amounts  of  GPS  data  a  mainframe  computer  should  be  used  rather  than  a 
microcomputer  to  reduce  processing  time. 
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APPENDIX  A.    FORTRAN  PROGRAM  GPSCON 


*  WRITTEN  BY  MA,  WEI-MING    05/08/88  * 

*  THIS  PROGRAM  READS  DATA  FROM  GEOID  DATA  FILE  AND  CALLS  * 

*  THE  SUBROUTINE  CONTOR  TO  PLOT  GEOID  HEIGHTS.  * 

*  THE  VARIABLES  USED  ARE:  * 


*  GEOHT 

*  NX 

*  NY 


GEOID  HEIGHTS  (m)  * 

NUMBER  OF  POINTS  IN  THE  X-DIRECTION  * 

NUMBER  OF  POINTS  IN  THE  Y-DIRECTION  * 

*  * 

REAL*4    GEOHT(22,15) 
DATA     NX,NY/22,15/ 
C 

C  READ  GEOID  HEIGHTS  FROM  GEOID  1  DATA  FILE 
C 

CALL  EXCMSOFILEDEF  8  DISK  GEOID  1  DATA  Al*) 
DO  200  J  =1,  NY 

DO  1001  =  NX,  1,-1 

READ(8,*)  GEOHT(IJ) 
100  CONTINUE 

200    CONTINUE 
C 

C  PLOT  THE  GEOID  HEIGHT  BY  CALLING  SUBROUTINE  CONTOR 
C 

CALL  CONTOR(GEOHT,NX,NY) 

STOP 

END 

SUBROUNTINE  CONTOR(A,NX,NY) 

*  * 

*  WRITTEN  BY  MA,  WEI-MING    05/08/88  * 

*  THIS  SUBROUTINE  CONTOURS  AN  NX  BY  NY  ARRAY  OF  * 

*  REGULARLY  SPACED  POINTS.  * 

*  NOTE:  THIS  ARRAY  MUST  BE  REAL*4.  * 

*  THE  VARIABLES  USED  ARE:  * 

*  A  :      SINGLE  PRECISION  NX  BY  NY  ARRAY  OF  * 

*  REGULARLYSPACED  POINTS.  * 


*  NX 

*  NY 

*  ZINC 


NUMBER  OF  POINTS  IN  THE  X-DIRECTION  * 

NUMBER  OF  POINTS  IN  THE  Y-DIRECTION  * 

CONTOUR  INTERVAL  * 

*  * 

DIMENSION  A(NX,NY) 

COMMON  WORK(5000) 
C 

C  SET  PARAMETERS  FOR  AXES 
C 
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XORIG  =  -55.0 

XSTP  =  2.0 

XMAX  =  -34.0 

YORIG  =  35.0 

YSTP  =  2.0 

YMAX  =  49.0 
C 

C  SET  CONTOUR  LEVEL 
C 

ZINC  =  0.1 

CALL  COMPRS 

CALL  SETCLR('CYAN') 
C 

C  SET  PAGE  AND  PLOT  SIZES,  SET  UP  AXES  FOR  PLOT 
C 

CALL  PAGE(8.0,6.0) 

CALLBCOMON(5000) 

CALL  AREA2D(7.0,5.0) 
C 

C  LABEL  AXES 
C 

CALL  XNAMECLONGITUDE  121  DEGREE  WEST  (MINUTES )$',1 00) 

CALL  YNAMECLATTTUDE  36  DEGREE  NORTH  (MINUTES )$',  100) 

CALL  GRAF(XORIG,XSTP,XMAX,YORIG,YSTP, YMAX) 

CALL  FRAME 
C 

C  MAKE  CONTOURS  AND  DRAW 
C 

CALL  SETCLR('RED') 

CALL  CONMIN(3.0) 

CALL  CONANG(60.) 

CALL  CONLIN(0,'MYCON','LABELS  *,  1,10) 

CALL  CONMAK(A,NX,NY,ZINC) 

CALL  CONTURO, 'LABELS'/DRAW) 
C 
C  END  PLOT 

CALLENDPL(O) 

CALLDONEPL 

RETURN 

END 
C  SUBROUTINE  MYCON(RARAY,IARAY) 

*  * 

*  THIS  SUBROUTINE  MAKES  NEGATIVE  CONTOURS  DASHED  AND  * 

*  THE  ZERO  LINE  HEAVYIER  * 

*  * 

DIMENSION  RARAY(2),IARAY(1) 
CALL  RESET(DASH') 
IF  (RARAY(l).GE.  0.)  GO  TO  10 
CALL  DASH 
10      RARAY(2)=1. 
IARAY(1)=1 
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IF(RARAY(l).EQ.O.)  IARAY(1)  =  2 

RETURN 

END 
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APPENDIX  B.    BASIC  PROGRAM  LOBS.BASIC 


5  REM  LEAST  SQUARES  BY  OBSERVATION 

10  DIML(N%,1) 

20  DIMA(N%,X%),  P(X%,X%),  EX(X%,1) 

30  DIM  CX(X%,1) 

40  DIM  AT(X%,N%),  R(50),  ATP(X%,X%) 

50  DIM  ATPA(X%,X%),  ATPL(X%,1),  AI(X%,X%) 

69  PRINT  "  INPUT  #  OF  OBSERVATIONS" 

70  NPUT  NO 

200  PRINT  "  INPUT  #  OF  PARAMETERS" 

270  INPUT  Nl 

350  ERASE  L 

351  ERASE  A 

352  ERASE  AT 

353  ERASE  P 

354  ERASE  ATP 

355  ERASE  ATPL 

356  ERASE  AI 

357  ERASE  CX 

358  ERASE  EX 

359  ERASE  ATPA 

380  DIM  L(NO,l),  A(NO,Nl),  AT(Nl,NO),  P(NO,NO),  ATP(Nl,NO), 

ATPL(N1,1),AI(N1,N1) 

390  DIM  CX(N1,1),  EX(NO,l),  ATPA(N1,N1) 

400  ITERO  =  0 

1 180  PRINT  "  INPUT  COEFFICIENTS  OF  OBSERVATION  EQUATIONS  AND 

WEIGHTS" 

1190  FORK=lTONO 

1200  FORJ  =  lTONl 

1220  PRINT  "COEFFICIENT  ",  K,  J 

1230  INPUT  A(K,J) 

1240  PRINT  A(K,J) 

1250  NEXT  J 

1260  PRINT  "  INPUT  OBSERVED  VALUE  ",  K 

1270  INPUT  L(K,1) 

1280  PRINT  "  INPUT  WEIGHT  OF  OBSERVED  VALUE  ",  K 

1300  INPUT   P(K,K) 

1390  NEXT  K 

1410  DOF  =  0 

1420  REM  LEAST  SQUARES 

1430  M  =  NO 

1450  L  =  Nl 

1460  PRINT  "  M  =",  M,  "Nl  =  ",  Nl,  "L  =  "L 

1470  FOR  I  =  1  TO  M 

1480  FORJ=lTOL 

1490  AT(J,I)  =  Aa,J) 

1500  NEXT  J 

1510  NEXT  I 

1520  REM  ATP  =  AT*P 

1530  N  =  NO 

1540  FORI=l  TOL 
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1550 

FORJ=l  TON 

1560 

ATP(IJ)  =  0 

1570 

FORK=l  TOM 

1580 

ATP(IJ)  =  ATP(I,J)  +  AT(I,J)  *  P(K,J) 

1590 

NEXTK 

1600 

NEXT  J 

1610 

NEXT  I 

1620 

REM  ATPA  =  ATP  *  A 

1630 

N  =  N1 

1640 

FOR  I  =  1  TO  L 

1650 

FOR  J  =  1  TO  N 

1660 

ATPA(I,J)  =  0 

1680 

ATPA(I.J)  =  0 

1690 

FOR  K  =  1  TO  M 

1710 

ATPA(IJ)  =  ATPA(I,J)  +  ATP(I,K)  *  A(K,J) 

1720 

NEXTK 

1730 

NEXT  J 

1740 

NEXT  I 

1750 

REM  ATPL  =  ATP  *  L 

1760 

N=  1 

1770 

FOR  I  =  1  TO  L 

1780 

FOR  J  =  1  TO  N 

1790 

ATPL(IJ)  =  0 

1800 

FORK=l  TOM 

1810 

ATPL(IJ)  =  ATPL(I,J)  +  ATP(I,K)  *  L(K,J) 

1820 

NEXTK 

1830 

NEXT  J 

1840 

NEXT  I 

1970 

REM  AI  =  INV(ATPA) 

1980 

I  =  N1 

1990 

M  =  N1 

2000 

N  =  I-  1 

2020 

MI  =  M-1 

2030 

FOR  J  =  1  TO  I 

2040 

FOR  K  =  1  TO  I 

2050 

AI(J,K)  =  ATPA(J,K) 

2055 

NEXTK 

2060 

NEXT  J 

2070 

FOR  K  =1  TO  I 

2080 

FOR  J  =  1  TO  MI 

2090 

R(J)  =  AI(1,J+1)/AI(1,1) 

2100 

NEXT  J 

2110 

R(M)  =  1!/AI(1,1) 

2120 

FOR  L  =  1  TO  N 

2130 

FOR  J  =  1  TO  MI 

2140 

AI(I,J)  =  AI(L+1,J+1)  -  AI(L+1,1)  *  R(J) 

2150 

NEXT  J 

2160 

AI(L,M)  =  -AI(L+1,1)  *  R(M) 

2170 

NEXTL 

2180 

FOR  J  =  1  TO  M 

2190 

AI(I,J)  =  R(J) 

2200 

NEXT  I 

2210 

NEXTK 
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2220 

REM  CX  =  AI  *  ATPL 

2230 

L  =  N1 

2240 

N=l 

2250 

M  =  N1 

2260 

FOR  I  =  1  TO  L 

2270 

FOR  J  =  1  TO  N 

2280 

CX(I,J)  =  O 

2300 

FOR  K  =  1  TO  M 

2310 

CX(I,J)  =  CX(I,J)  +  AI(I,K) 

*  ATPL  (K,J) 

2320 

NEXTK 

2325 

PRINT  "  PARAMETER  # 

",  I,  CX(I,J) 

2330 

NEXT  J 

2340 

NEXT  I 

2480 

PRINT  "  RESIDUAL" 

2490 

L  =  NO 

2500 

N=l 

2510 

M  =  N1 

2520 

REM  EX  =  A  *  CX 

2530 

FOR  I  =  1  TO  L 

2540 

FOR  J  =  1  TO  N 

2560 

EX(IJ)  =  0 

2570 

FOR  K  =  1  TO  M 

2580 

EX(IJ)  =  EX(I,J)  +  A(I,K) 

*  CX(KJ) 

2590 

NEXTK 

2600 

NEXT  J 

2610 

NEXT  I 

2620 

SD  =  0 

2630 

FOR  K  =  1  TO  NO 

2640 

V  =  EX(K,1)  -  L(K,1) 

2650 

PRINT  K,  V 

2660 

SD  =  SD  +  V  *  V 

2670 

NEXTK 

2680 

SD  =  SQR(SD  /  (NO  -  Nl  + 

DOF)) 

2690 

PRINT  SD 

2695 

PRINT  "  STD.DEV",  SD 

2700 

ITER  =  ITER  +  1 

2710 

IF  ITER  <  3  GOTO  350 

2720 

END 
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APPENDIX  C.    FORTRAN  PROGRAM  GPSDIS 


*  WRITTEN  BY  WEI-MING  MA  05/14/88  * 

*  THIS  PROGRAM  READS  DATA  FROM  TERMINAL  AND  COMPUTES  * 

*  THE  ORTHOMETRIC  HEIGHTS  FOR  THE  CHECK  MARKS  THEN  * 

*  COMPUTES  THE  DIFFERENCE.  * 

*  THE  VARIABLES  USED  ARE:  * 

*  RX,RY,RZ  :      THE  COORDINATES  OF  THE  REFERENCE  MARK  (m)  * 

THE  COORDINATES  OF  THE  CHECK  MARK  (m)  * 

THE  COORDINATES  DIFFERENCE  (m)  * 

THE  ORTHOMETRIC  HEIGHT  OF  CHECK  MARK  (m)  * 

THE  LEVELED  ORTHOMETRIC  HEIGHT  (m)  * 

THE  ELLIPSOID  HEIGHT  DIFFERENCE  (m)  * 

THE  GLOBAL  GEOID  HEIGHT  OF  CHECK  MARK  (m)  * 

THE  LEVELED  ORTHOMETRIC  HEIGHT  (m)  * 

THE  DIFFERENCE  BETWEEN  H  AND  LH  (cm)  * 

THE  COEFFICIENTS  OF  5-PARAMETER  GEOID  MODEL  * 

THE  COEFFICIENTS  OF  6-PARAMETER  GEOID  MODEL  * 


*  CX,CY,CZ 

*  DX,DY,DZ 

*  H 

*  LH 

*  DH 

*  GN 

*  LH 

*  DIF 


* 


*  A,B,C,D 

*  A1,B1,C1 

*  D1,E1 

*  IPAJPA1    :      THE  NUMBER  OF  PARAMETERS  * 

*  IW          :      THE  OUTPUT  UNIT  * 

*********************************************************************** 

REAL*8  RX,  RY,  RZ,  CX,  CY,  CZ,  DX,  DY,  DZ,  LH,  GN, 

REAL*8  A,  B,  C,  D,  Al,  Bl,  CI,  Dl,  El,  H,  DH,  DIF, 

INTEGER  IPA,  IPA1,IW 
CHARACTER*  1  RESPN1,  RESPN2,  RESPN3 
IW  =  6 
IPA  =  0 
IPA1=0 
C 

C  READ  DATA  FROM  TERMINAL 
C 
10      PRINT  *,  'ENTER  #  OF  PARAMETERS' 
READ  *  IPA 
IF  (IPA  .EQ.  5  .OR.  IPA  .EQ.  6)  THEN 

PRINT  *,  'ENTER  THE  X  COORDINATE  IN  WGS  84  OF  THE  ', 
&  'REFERENCE  MARK' 

READ  *,  RX 

PRINT  *,  'ENTER  THE  Y  COORDINATE  IN  WGS  84  OF  THE  ', 
&  'REFERENCE  MARK' 

READ  *,  RY 

PRINT  *,  'ENTER  THE  Z  COORDINATE  IN  WGS  84  OF  THE  ', 
&  'REFERENCE  MARK' 

READ  *,  RZ 
ELSE 

STOP 
END  IF 
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IF  (RESPN1  .EQ.  'Y')  THEN 

PRINT  *,*DO  YOU  WANT  TO  CHANGE  THE  CHECH  MARK  ("Y"  OR  ', 
&  "'N'*)' 

READ  *,  RESPN3 
IF  (RESPN3  .EQ.  'Y')  GO  TO  30 
GO  TO  40 
END  IF 
20      PRINT  *,'DO  YOU  WANT  TO  CHANGE  THE  REFERENCE  MARKS  ("Y" ', 
&  'OR    '¥')' 

READ  *,  RESPN1 
IF  (RESPN1  .EQ.  *Y')  GO  TO  10 
30      PRINT  VENTER  THE  X  COORDINATE  IN  WGS  84  OF  THE  CHECK.', 
&  'MARK' 

READ  *,  CX 

PRINT  *,'ENTER  THE  Y  COORDINATE  IN  WGS  84  OF  THE  CHECK  MARK' 
READ  *  CY 

PRINT  VENTER  THE  Z  COORDINATE  IN  WGS  84  OF  THE  CHECK  MARK' 
READ  *,  CZ 

PRINT  VENTER  THE  ELLIPSOID  DIFFERENCE,  DH' 
READ  *,  DH 

PRINT  VENTER  THE  GLOBAL  GEOID  HEIGHT  OF  THE  CHECK  MARK', 
&  ',  GN' 

READ  *,  GN 

PRINT  VENTER  THE  LEVELED  OTHOMETRIC  HEIGHT  OF  THE  CHECK  ' 
&  ,  'MARK' 

READ  *,  LH 
DX  =  CX  -  RX 
DY  =  CY  -  RY 
DZ  =  CZ-  RZ 

IF  (IPA  .EQ.  IPA1  .AND.  RESPN1  .EQ.  *N'  .AND.  IPA  .EQ.  5)  GO  TO  50 
IF  (IPA  .EQ.  IPA1  .AND.  RESPN1  .EQ.  'N'  .AND.  IPA  .EQ.  6)  GO  TO  60 
40      IF  (IPA  .EQ.  5)  THEN 

PRINT  *, 'ENTER  THE  ELLIPSOID  HEIGHT  HO' 
READ  *,  HO 

PRINT  VENTER  THE  COEFFICIENT  A' 
READ*,  A 

PRINT  VENTER  THE  COEFFICIENT  B' 
READ  *,  B 

PRINT  VENTER  THE  COEFFICIENT  C 
READ  *,  C 

PRINT  VENTER  THE  COEFFICIENT  D' 
READ*,D 
C 

C  THE  5-PARAMETER  GEOID  MODEL 
C 
50  H  =  -GN  +  DH  +  HO  +  A*DY  +  B*DX**2  +  C*DY**2  +  D*DX*DY 

ELSE 

PRINT  VENTER  THE  ELLIPSOID  HEIGHT  HO' 

READ  *,  HO 

PRINT  VENTER  THE  COEFFICIENT  A*" 

READ  *,  Al 

PRINT  VENTER  THE  COEFFICIENT  B"' 

READ*,B1 
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PRINT  *, 'ENTER  THE  COEFFICIENT  C" 
READ  *  CI 

PRINT  *,'ENTER  THE  COEFFICIENT  D*" 
READ  *,  Dl 

PRINT  *,'ENTER  THE  COEFFICIENT  Em 
READ  *,  El 
C 

C  THE  6-PARAMETER  GEOID  MODEL 
C 
60  H  =  -GN  +  DH  +  HO  +  A1*DX  +  B1*DZ  +  C1*DX**2  +  D1*DY**2 

&  +  E1*DX*DY 

END  IF 
C 

C  COMPUTE  THE  DIFFERENCE  AND  PRINT  THE  RESULT 
C 

DIF  =  (  H  -  LH)  *  100. 
WPJTE(IW,  1)IPA 
WRITE(IW,  2)  LH,  H,  DIF 

1  FORMAT  (5X,  II, '  PARAMETERS') 

2  FORMAT  (/5X,'H-LEVELBSfG    =',  F10.3,  '  M* 
&  /5X,'H-GPS  =',  F10.3,  *  M' 

&  /5X,'DIFFERENCE    =',  F10.3,  '  CM'/) 

IPA1=IPA 
H  =  0. 
DIF  =  0. 

PRINT  *,  'MORE  COMPUTATION  ("Y"  OR  "N")' 
READ  *,  RESPN2 
IF  (RESPN2  .EQ.  'Y' )  GO  TO  20 
STOP 
END 
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APPENDIX  D.    FORTRAN  PROGRAM  DISTCO 


*  WRITTEN  BY  MA,  WEI_MING    05/11/88  * 

*  THIS  PROGRAM  COMPUTES  THE  DISTANCE  BETWEEN  THE  BENCH      * 

*  MARKS  THEN  CALLS  THE  SUBROUTINE  C  TO  CALCULATE  THE  * 

*  COVARIANCE  BETWEEN  THE  BENCH  MARKS.  * 

*  THE  VARIABLES  USED  ARE:  * 

*  X,Y,Z      :      THE  COORDINATES  OF  X,Y,Z,  IN  WGS  84  (m)  * 

*  DIS         :      THE  DISTANCES  BETWEEN  THE  MARKS  (m)  * 

*  CQ         :      THE  CO  VARIANCES  BETWEEN  THE  BENCH  MARKS     * 

*  * 

REAL*8  X(10),  Y(10),  Z(10) 

REAL*8  DIS(10,10),  CQ(10,10) 

DATA     DIS/100*0./,  CQ/100*0./ 

IW  =  9 

N  =  7 
C 

C  READ  X,  Y,  Z  FROM  THE  DATA  FILE 
C 

CALL  EXCMSC'FILEDEF  8  DISK  MBXYZ  DATA  Al') 

DO100I  =  l,N 

READ(8,*)  X(I),  Y(I),  Z(I)) 
100    CONTINUE 
C 

C  COMPUTE  THE  DISTANCE  AND  PRINT  THE  RESULTS 
C 

DO300I=l,N 
WRITE(IW,1)I 

1  FORMAT(/lX,'FROM  BENCH  MARK  #',  12,'    TO',  'DISTANCE  (m)7 
&  IX,  53('-')) 

DO200J  =  l,N 

IF  (J  .EQ.  I)  GO  TO  10 

DIS(IJ)  =  SQRT((X(J)  -  X(I))  **  2  +  (Y(J)  -  Y(I))  **  2 
&  +(Z(J)-Z(I))**2) 

10  WRITE(IW,2)  J,  DIS(I,J) 

2  FORMAT(llX,  12,  14X,  G20.8) 
200          CONTINUE 

300    CONTINUE 
C 

C  CALL  SUBROUTINE  C  TO  CALCULATE  THE  CQ  AND  PRINT  THE  RESULTS 
C 

WRITE(IW,3) 

3  FORMAT(/3X,  'CQ  =  ') 
CALL  C(DIS,N,CQ) 
DO400J=l,N 

WRITE(IW,4)  (CQ(I,J),  I  =  1,  N) 

4  FORMAT(3X,  7(G16.7,  2X) 
400    CONTINUE 
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STOP 
END 

SUBROUNTINE  C(DIS,N,CQ) 

********************************************************************** 

*  * 

*  WRITTEN  BY  MA,  WEI_MING    05/11/88  * 

*  THIS  SUBROUTINE  COMPUTES  THE  COVARIANCES  OF  THE  RANDOM  * 

*  ERROR  DISTRIBUTION  FUNCTION:  * 

*  C(R)  =  A  +  B  *  SIN(D*DIS)  * 

*  THE  VARIABLES  USED  ARE:  * 


*  A 

*  B 

*  C 

*  DIS 


STANDARD  DEVIATION  OF  THE  OBSERVATIONS  * 

CONSTANT  * 

CONSTANT  (RADIANS/METERS)  * 

DISTANCE  (METERS)  * 


REAL*8  DIS(10,10),  CQ(10,10),  A,  B,  D 

A  =  0.137153 

B  =  0.147951 

D  =  2.85958 

DO200J  =  l,N 
DO100I=l,N 

CQ(IJ)  =  A  -  B  *  SIN(D  *  DIS(IJ)  /  33300.) 
100         CONTINUE 
200    CONTINUE 

RETURN 

END 


76 


REFERENCES 


Ashjaee,  J.,  New  results  on  the  accuracy  of  the  C/A  code  GPS  receiver.  In: 

Proceedings  of  the  First  International  Symposium  on  Precise  Positioning  with 

the  Global  Positioning  System,  vol.  1,  pp.  207-214,  U.S.  Department  of 

Commerce,  Rockville,  Maryland,  April  15  to  19,  1985. 
Baker,  P.  J.,  Global  Positioning  System  (GPS)  policy.  In:  Proceedings  of  the  Fourth 

International  Geodetic  Symposium  on  Satellite  Positioning,  vol.  1,  pp.  51-64, 

Sponsored  by  the  Defense  Mapping  Agency  and  the  National  Geodetic  Survey, 

at  Austin,  Texas,  April  28  to  May  2, 1986. 
Bodnar,  A.  N.  JR,  User's  Guide  for  the  Establishment  of  Tidal  Bench  Marks  And 

Levelling  Requirements  for  Tide  Stations,  p.  1 1 ,  National  Geodetic  Survey 

Charting  and  Geodetic  Services  National  Ocean  Service,  NOAA,  Rockville, 

Maryland,  December  8,  1975. 
Bomford,  G.,  Geodesy.  4th  ed.,  pp.  736-739,  Clarendon  press,  Oxford,  1980. 
Bouchard,  R.  H.,  Optimized  Observation  Periods  Required  to  Achieve  Geodetic 

Accuracies  Using  the  Global  Positioning  System.  M.  S.  Thesis,  Naval 

Postgraduate  School,  Monterey,  California,  March  1988. 
Davis,  R.  E.,  F.  S.  Foote,  J.  M.  Anderson,  and  E.  M.  Mikhail,  Surveying:  Theory 

and  Practice.  6th  ed.,  pp.  160-168,  McGraw-Hill,  New  York,  1986. 
Defense  Mapping  Agency,  Department  of  Defense  World  Geodetic  System  1984, 

DMA  Technical  Report  8350.2,  pp.  3-10,  3-1 1,  Washington,  D.C.,  1987. 
Defense  Mapping  Agency,  Geodesy  for  the  Layman.  DMA  Technical  Report,  80- 

003,  pp.  24,  39,  64,  December  1983. 
Denker,  H.  and  Wenzel,  G.,  Local  Geoid  Determination  and  Comparison  with  GPS 

Result,  Bulletin  G6od6sique.  61(4),  349-366,  pp.  349-366,  1987. 
Ewing,  C.  E.  and  M.  M.  Mitchell,  Introduction  to  Geodesy  pp.  9,  176-180,  Elsevier 

North-Holland,  Inc.  1976. 
Fell,  P.  J.,  The  Use  of  Standard  Values  and  Refraction  Bias  Parameters  in  Orbit 

Determination,  The  Canadian  Surveyor.  29(3),  pp.  301-305,  September  1975. 
Jeyapalan,  K.,  Calibration  of  comparators  by  the  method  of  collocation,  pp.  2-7, 

Unpublished  Report.  Topographic  Division,  U.S.  Geological  Survey,  Reston, 

Virginia,  August  1977. 


77 


Kaula,  W.  M.,  The  Need  for  Vertical  Control,  Surveying  and  Mapping.  47(1),  pp. 

57-64,  March  1987. 
King,  R.  W.,  E.  G.  Masters,  C.  Rizos,  and  J.  Collins,  Surveying  with  GPS,  p.  128, 

School  of  Surveying,  The  University  of  New  South  Wales,  Kensington  N.S.W, 

Australia,  1985. 
Mikhail,  E.  M.,  Observations  and  Least  Squares,  pp.  418-426,  IEP-Dun-Donnelley, 

New  York,  1976. 
Milliken,  R.  J.,  and  C.  J.  Zoller,  Principle  of  operation  of  NAVSTAR  and  System 

Characteristics.  In:  Global  Positioning  System,  vol.  1,  pp.  3-14,  Institute  of 

Navigation,  Washington,  D.C.,  1980. 
Moritz,  H.  and  W.  A.  Heiskanen,  Physical  Geodesy,  pp.  82-86,  Reprint,  Institute  of 

Physical  Geodesy,  Technical  University,  Graz,  Austria,  1984. 
NASA,  Directory  of  Station  Locations.  4th  ed.,  p.  1-11,  Godard  Space  Flight 

Center,  Greenbelt,  Maryland,  February  1978. 
Rapp,  R.  H.,  and  J.  Y.  Cruz,  Spherical  Harmonic  Expansions  of  the  Earth's 

Gravitational  Potential  to  Degree  360  Using  30'  Mean  Anomalies,  Reports  of 

the  Department  of  Geodetic  Science  and  Surveying.  The  Ohio  State  University, 

Columbus,  Ohio,  Report  No.376,  pp.  1-2,  December  1986. 
Remondi,  B.  W.,  Global  Positioning  System  carrier  phase:  description,  and  use, 

Bulletin  Gdod6sique.  59(4),  361-377, 1985. 
Remondi,  B.  W.,  Using  the  Global  Positioning  System  (GPS)  Phase  Observable  for 

Relative  Geodesy:  Modelling.  Processing  and  Results.  Ph.  D.  Dissertation,  The 

University  of  Texas  at  Austin,  May  1984. 
Torge  W.,  Geodesy,  pp.  46,  135,  138,  Walter  de  Gruyter,  New  York,  1980. 
Trimble  Navigation,  Trimble  Model  4000SX    GPS    Surveyor-Preliminary- 

Installation  and  Operation  Manual  (Revision.  8/1/87),  pp.  96,  Sunnyvale, 

California,  1987a. 
Trimble  Navigation,  Trimvec  GPS  Survey  Software  Preliminary  User's  Manual. 

Revision  B,  pp.  62,  Sunnyvale,  California.,  1987b. 
Tziavos,  I.  N.,  Determination  of  Geoidal  heights  and  deflecting  of  vertical  for  the 

Hellenic  area  using  Heterogeneous  data,  Bulletin  Gdoddsique.  61(2),  pp.  177- 

197,  1987. 
U.S.  Department  of  Commerce,  Coast  and  Geodetic  Survey,  Washington,  D.C., 

Vertical  Control  Data.  Sea-Level  Datum  of  1929,  November  1961. 


78 


Wells,  D.  E.,  Guide  to  GPS  Positioning.  Canadian  GPS  Associates,  Fredericton, 

New  Brunswick,  Canada,  1986. 
Wolf,  P.  P.,  Adjustment  Computations:  Practical  Least  Squares  for  Surveyors.  2nd 

ed.,  3rd  printing,  p.  91,  Landmark  Enterprises,  California,  1980. 
Zilkoski,  D.  B.,  GPS  Satellite  Surveys  and  Vertical  Control,  Paper  presented  at  the 

GPS-88  Engineering  Applications  of  GPS  Satellite  Surveying  Technology  at 

Nashville,  Tennessee,  May  1 1-14,  1988. 


79 


INITIAL  DISTRIBUTION  LIST 


No.  Copies 


1 .  Defense  Technical  Information  Center  2 
Cameron  Station 

Alexandria,  VA  22304-6145 

2.  Library,  Code  0142  2 
Naval  Postgraduate  School 

Monterey,  CA  93943-5002 

3.  Chairman  (Code  68Co)  1 
Department  of  Oceanography 

Naval  Postgraduate  School 
Monterey,  CA  93943 

4.  Director  (Code  HO)  1 
Defense  Mapping  Agency  Hydrographic 

Topographic  Center 
6500  Brookes  Lane 
Washington,  D.C.  20315 

5.  Director,  Charting  and  Geodetic  1 
Services  (N/CG) 

National  Ocean  and  Atmospheric 

Administration 

Rockville,  MD  20852 

6.  Prof.  Stevens  P.  Tucker  2 
Department  of  Oceanography  (Code  68Tx) 

Naval  Postgraduate  School 
Monterey,  CA  93943 

7.  Prof.  Kandiah  Jeyapalan  2 
Department  of  Civil  Engineering 

Iowa  State  University 
Ames,IA  50011 


80 


8.  Library  of  Chinese  Naval  Academy 
P.O.  Box  8494  Tso-Ying, 
Kaohsiung,  Taiwan 

Republic  of  China 

9.  Library  of  Chung-Cheng  Institute  of  Technology 
Tashih,  Tao-Yuan,  Taiwan 

Republic  of  China 

10.  Chinese  Naval  Hydrographic  and 
Oceanographic  office 

P.O.  Box  8505  Tso-Ying, 
Kaohsiung,  Taiwan 
Republic  of  China 

1 1 .  Director,  Chinese  Naval  Hydrographic  and 
Oceanographic  Office 

P.O.  Box  8505  Tso-Ying, 
Kaohsiung,  Taiwan 
Republic  of  China 

12.  Chairman,  Survey  Engineering  Department 
Chung-Cheng  Institute  of  Technology 
Tashih,  Tao-Yuan,  Taiwan 

Republic  of  China 

13.  Prof.  Von  Schwind 

Department  of  Oceanography  (Code  68Vs) 
Naval  Postgraduate  School 
Monterey,  CA  93943 

14.  CDR.  Kurt.  J.  Schnebele 

Department  of  Oceanography  (Code  68Sn) 
Naval  Postgraduate  School 
Monterey,  CA  93943 

15.  LT.  Richard  H.  Bouchard 
Navoceancomcen  /  JTWC 
Comnavmarianas  Box  12 

FPO  San  Francisco  96630-2926 


81 


16.  Yu,Ta-Te 

Chinese  Naval  Hydrographic  and 
Oceanographic  Office 
P.O.  Box  8505  Tso-Ying, 
Kaohsiung,  Taiwan 
Republic  of  China 

17.  Chang,  Chin-Wen 

Chinese  Naval  Hydrographic  and 
Oceanographic  Office 
P.O.  Box  8505  Tso-Ying, 
Kaohsiung,  Taiwan 
Republic  of  China 

18.  Wang,  Chih-Ping 

Chinese  Naval  Hydrographic  and 
Oceanographic  Office 
P.O.  Box  8505  Tso-Ying, 
Kaohsiung,  Taiwan 
Republic  of  China 

1 9.  The  Hydrographer  of  Pakistan 
Hydrographic  Department 
Naval  Headquarters 
Islamabad,  Pakistan 

20.  Dr.  Muneendra  Kumar 
Defense  Mapping  Agency 
Hydrographic/Topographic  Center 
6500  Brooks  Lane 
Washington,  D.C.  20315-0030 

21.  Dr.  Narendra  K.  Saxena 
Department  of  Civil  Engineering 
University  of  Hawaii  at  Manoa 
2540  Dole  Street,  Holmes  383 
Honolulu,  HI  96822 

22.  Ma,  Wei-Ming 

Fl.  4  No.  5  Ln.  1558  Yu-Chen  Rd., 
Kaohsiung,  Taiwan  80402 
Republic  of  China 


82 


23.    Mr.  Gary  Fredrick 
NOAA,  NOS 
Code  N/MOP  222 
Bin  C15700 

7600  Sand  Point  Way,  NE 
Seattle,  WA  98115-0070 


83 


Local  geoid  determina- 
tion using  the  Global 
Positioning  System. 


II   JAX  93 


38287 


Thpsis 

MICH   Ma 

c.l      Local  geoid  determina- 
tion using  the  Global 
Positioning  System. 


